#include <string.h>

#include "investment_include.h"

int main(void)
 {
	 num=100;
	 avg_LS_MC=0;
	 contQI=0;
	 contalfasup=0;
	 contQtot1=0;
	 contINN=0;
	 contIMM=0;
	 contfront=0;
	 contfront2=0;
	 contA1entry=0;
	 contA1front=0;
	 contFE=0;
	 DebonGDPm1=0;
	 DebonGDPm2=0;
	 contSI=0;
	 contMONOP=0;
	 contcpi1=0;
	 contppi1=0;
	 contCC=0;
	 contIMMnaz=0;
	 contIMMfor=0;
	 umin=0;
			
	 for (imc=1; imc<=MC; imc++)
	 {
		 cout << "\n\n\n --------------- C I C L O   MC   N:   " << num << " ---------------" << endl;   
		 if (imc<=3)
		 {
			 GENFILEOUTPUT("output",num);
			 GENFILEOUTPUT1("out",num);
			 GENFILEPROD1("A1_",num);
			 GENFILEPROD2("A2_",num);
			 GENFILEPRODALL1("A1all_",num);
			 GENFILEPRODALL2("A2all_",num);
		 }
		 		
		 GENFILEYMC("ymc",num);
		 INITIALIZE();						// Fare check seed!
		 			 
		 for (t=1; t<=T; t++)   
		 {
			 if (echo==1)
				 cout << "\n\n\n T E M P O    T = " << t << endl;
			 
			 if (imc<=3 && t >= 500)
			 {
					 nsize=imc*1000+t-500+1;
					 GENFILESIZE1("ZS1_",nsize);
					 GENFILESIZE2("ZS2_",nsize);
			 }
					 
			 //cout << "\nqua1" << endl;
			 MACH();
			 //cout << "\nqua2" << endl;
		 	 BROCHURE();
			 //cout << "\nqua3" << endl;
		 	 INVEST();		// calcolo investimenti
			 //cout << "\nqua4" << endl;
		 	 PRODMACH();
			 //cout << "\nqua5" << endl;
		 	 COMPET2();		// calcolo competitivit‚àö‚Ä†, market share 
			 //cout << "\nqua6" << endl;
		 	 PROFIT();		// calcolo domanda di lavoro, profitti e ricchezza netta imprese
			 //cout << "\nqua7" << endl;
		 	 MACRO();		// calcolo variabili aggregate
			 //cout << "\nqua8" << endl;
		 	 ENTRYEXIT();	// Entry & Exit
			 //cout << "\nqtua9" << endl;
			 if (flagTCGEN==1)
				 TECHANGEX();	// Technical change exogenous
			 else
			 {
				 if (flagTCGEN==2)
					 TECHANGEND();  // Technical change endogenous
				 else cerr << "\n ERRORE flagTCGEN!" << endl;
			 }
			 //cout << "\nqua10" << endl;
		 	 SAVE();		// Salva risultati su file
			 //cout << "\nqua11" << endl;
		 	 UPDATE();		// aggiorna variabili per prossimo periodo
		 	 OVERBOOST();
		 	 
			 if (t==1) 
			 {
			 	 Y0=Y;
			 	 Am0=Am(1);
				 LD0=LD;
				 Pi20=Pitot2;
				 Pi10=Pitot1;
			 }	
			 if (parbertot > parbermax)
			 {
				 parbermax=parbertot;
				 tbermax=t;
			 }
			 if (parbertot1 > parbermax1)
			 {
				 parbermax1=parbertot1;
				 tbermax1=t;
			 }
			 
		} 
	 	 debug.close();
		
		 if (flagbug == 0)
		 {
		 	 avg_prod=(log(Am(1))-log(Am0))/(T);
		 	 avg_y=(log(Y)-log(Y0))/(T);
			 avg_LD=(log(LD)-log(LD0))/(T);
			 avg_Pi2=(log(Pitot2)-log(Pi20))/(T);
             avg_Pi1=(log(Pitot1)-log(Pi10))/(T);
			 Umean/=(T);
			 
		 	 cout << "\n\n Avg. prod. (100-600)  :   " << avg_prod << endl;
		 	 cout << "\n Avg. Y     (100-600):   " << avg_y << endl;
		 	 cout << "\n Avg. LD     (100-600):   " << avg_LD << endl;
			 if (flagEXTRA==1)
			 {
				 avg_LS=(log(LS)-log(LS0))/T;
				 avg_LS_MC+=avg_LS;
				 cout << "\n Avg. LS         :     " << avg_LS << endl;
			 }
			 cout << "\n Disoccupazione: (1-600): " << Umean << endl;
			 cout << "\n Avg. Pi2     (100-600):   " << avg_Pi2 << endl;
			 cout << " Avg. Pi1     (100-600):   " << avg_Pi1 << endl;
			 cout << "\n N. max. uscite industria 1:   "<< nextmax1 << endl;
		 	 cout << " N. max uscite industria 2:   " << nextmax2 << endl;
			 
			 if (flagTCGEN==2)
			 {
				 cout << "\n Quota max imprese entranti industria 2: " << f2max << endl;
				 cout << "\n Max par. Bernulli INN: " << parbermax << endl;
				 cout << " Tempo par. Ber. max INN: " << tbermax << endl;
				 cout << " Max Par. Bernulli IMM: " << parbermax1 << endl;
				 cout << " Tempo par. Ber. max IMM: " << tbermax1 << endl;
				 if (flagFRONT==1)
				 {
					 if (flagFRONTEX < 3)
					 {
						 cout << "\n Produttivita Frontiera : " << A1front(1) << endl;
						 cout << " Produttiivta max incumb.: " << A1max << endl;
						 if (A1max==A1front(1))
							 contA1front++;
						 cout << " Produttivita Frontiera pr. : " << A1pfront(1) << endl;
						 cout << " Produttiivta max incumb. pr.: " << A1pmax << endl;
						 if (A1pmax==A1pfront(1))
							 contA1front++;
					 }
					 else
					 {
						 cout << "\n Produttivita Max Frontiera : " << A1fmax << endl;
						 cout << " Produttivita Min Frontiera : " << A1fmin << endl;
						 cout << " Produttività Media Frontiera : " << A1favg << endl;
						 cout << " Deviazione Standard Produttività Frontiera : " << A1fstd << endl;
						 cout << " Produttiivta max incumb.: " << A1max << endl;
						 if (A1max>=A1fmax)
							 contA1front++;
						 cout << " Produttivita Max Frontiera P : " << A1pfmax << endl;
						 cout << " Produttivita Min Frontiera P : " << A1pfmin << endl;
						 cout << " Produttività Media Frontiera : " << A1pfavg << endl;
						 cout << " Deviazione Standard Produttività Frontiera : " << A1pfstd << endl;
						 cout << " Produttiivta max incumb. P: " << A1pmax << endl;
						 if (A1pmax>=A1pfmax)
							 contA1front++;
						 //A1ratio=Am1bis/A1favg;
						 //A1pratio=Am1/A1pfavg;
						 A1ratio=A1avg/A1favg;
						 A1pratio=A1pavg/A1pfavg;
						 cout << " Prod. developing / Prod. developed : " << A1ratio << endl;
						 cout << " Prod. P developing / Prod. P developed : " << A1pratio << endl;
					 }
				 }
				 cout << " Produttivita totale: " << Am(1) << endl;
				 cout << "\n Wtot1: " << Wtot1 << endl;
				 cout << " Wtot2: " << Wtot2 << endl;
			 }
			 
			 if (aliq>0 || flagC>0)
			 {
				 cout << "\n def/GDP: " << DefonGDP;
				 cout << "\n deb/GDP: " << DebonGDP;
				 DebonGDPm1+=DebonGDP;
				 DebonGDPm2+=DebonGDP*DebonGDP;
			 }			 			 
		 }		 		 
	 	 num++;		 
		 cout << "\n\n freq, imprese che investono ma non riescono a produrre: " << contQI << endl;
		 cout << " freq. alfa > alfasup: " << contalfasup << endl;
		 cout << " freq. I=0:" << contQtot1 << endl;
		 cout << " freq. INN: " << contINN << endl;
		 cout << " freq. IMM: " << contIMM << endl;
		 cout << " freq. IMM nazionali: " << contIMMnaz << endl;
		 cout << " freq. IMM foreign: " << contIMMfor << endl;
		 cout << " freq. A1 > A1front: " << contfront << endl;
		 cout << " freq. A1p > A1pfront: " << contfront2 << endl;
		 cout << " freq. A1entry > A1 incumbent: " << contA1entry << endl;
		 cout << " freq. piena occupazione: " << contFE << endl;
		 cout << " freq. impresa leader settore 1 ha quota > .99: " << contMONOP << endl;
		 cout << " freq. Cpast/Creal > .50: " << contCC << endl;
		 cout << " freq. p1 < " << pmin << ": " << contppi1 << endl;
		 cout << " freq. p2 < " << pmin << ": " << contcpi1 << endl;
		 cout << " disoccupazione minima (0 default): " << umin << endl;
		 if (flagEXTRA==1)
			 cout << " \nTasso di crescita popolazione MC: " << avg_LS_MC/MC << endl; 
	 }
	 DebonGDPm1/=MC;
	 DebonGDPm2=sqrt((DebonGDPm2/MC)-(DebonGDPm1*DebonGDPm1));
	 cout << "\n\n MC - Deb/GDP - mean: " << DebonGDPm1 << "\t\t std. dev.: " << DebonGDPm2 << endl;
	 cout << "\n n. repliazioni MC che terminano con imprese su frontiera: " << contA1front << endl;
	 cout << "\n Percentuali di turni in cui l'economia è in pieno impiego: " << contFE/T/MC << "\n\n" << endl;

	 return(0);
}

//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------

// FUNCTION BODIES

//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------

void INITIALIZE(void)
{
	//cout << "\nInitialization..." << endl;
	
	ofstream inv_ymc(nomefile2);
	
	if (imc == 1)
	{
		seed=-2;							// Inizializzo seed per numeri casuali normale
		p_seed=&seed;						// Puntatore
//		cout << "puntatore random: " << p_seed;
		srand(seed);
	}
	
	if (imc <=3)
		INTFILE();	// Scrivo intestazione dei file

	
	if (alfa > alfasup)
		cerr << "\nERRORE: alfa > alfasup!" << endl;
	
	Ir2=1; 
	w=w0;								// Inizializzo salario
	wstar=w0;
	Am(2)=A0;							// Inizializzo produttivit‚àö‚Ä† media
	A1=A0;								// Inizializzo produttivit‚àö‚Ä† imprese macchinari
	A1p=A0;
	A2=A0;
	
	cpi(2)=(1+mi2)*w0/A0;				// Inizializzo cpi
	ppi(2)=(1+mi1)*w0/A0;
	A=A0;								// Inizializzo produttivit‚àö‚Ä† macchinari
	C=w0/A0;							// Inizializzo costi macchine: C=w/A
	c2=w0/A0;
	c1=w0/(A0*a);
	p2=(1+mi2)*w0/A0;
	p1=(1+mi1)*w0/(A0*a);
	
	Em2=1;								// delle 2 industrie
	l2=1;								// Uguali ad 1 per via del logaritmo
	W1=W10;
	W2=W20;
	K=K0;								// Inizializzo capitale impresa
	LS=LS0;
	U(2)=1;
	N1r=double(N1);
	N2r=double(N2);
	E2=1;
	f1=1/N1r;
	f2=1/N2r;
	ueff=0;
	next1=0;
	next2=0;
	t0=1;
	
	if (flagC==1)
		G=0.005*LS0;
	else G=0;
	
	Deb1=0;
	Deb2=0;
	nextmax1=0;
	nextmax2=0;
	diff_w=0;
	diff_cpi=0;
	i=0;
	EI=0;
	SI=dim_mach;						// Scrappo una macchina per turno
	I=dim_mach;
	
	if (flagTC==2 && flagTCGEN==1)		// Con progresso tecnico deterministico inizializzo
	{									// per scrappare tutte macchine ogni turno
		Am=1/1.005;
		SI=K0;
		I=K0;
		A1=A0*(1+uu);
	}
	
	Q1=0;
	mat2d=0;
	Q2=0;
	Divtot=0;
	die=0;
	Cpast=0;
	mu2=mi2;
	Netot=0;
	flagA2=0;
	flagA1=0;
	Gtot=0;
	GonGDP=0;
	parbermax=0;
	parbermax1=0;
	Ke=0;
	Kexp=0;
	f2max=0;
	Td[0]=0;
	Pat=0;
	Patw=0;
	Umean=0;
	size1=0;
	size2=0;
	//if (aliq>0)
	//	Deb=-w0*LS0*wu*0.9567;
	//else
		Deb=0;
	
	step=N2/N1;							// Inizializzo matrice matching macchine
	cont=0;
	
	if (N2%N1 != 0)
		cerr << "\n\n ERRORE DI INIZIALIZZAZIONE!!!!!" << endl;
	Match=0;
	
	for (i=1; i<=N1; i++)
	{
		cont+=step;
		for (j=0; j<step; j++)
		{
			Match(cont-j,i)=1;
			fornit(cont-j)=i;
		}
	}
	
	for (tt=1; tt<=T; tt++)				
	{ 
		for (i=1; i<=N1; i++)   
		{ 
			for (j=1; j<=N2; j++)
			{
				g[tt-1][i-1][j-1]=0; 
				gtemp[tt-1][i-1][j-1]=0;
				g_c[tt-1][i-1][j-1]=0;
				age[tt-1][i-1][j-1]=0;
			}
		}
	}
	
	for (j=1; j<=N2; j++)					// Inizializzo matrice 3d frequenze macchine
	{ 	 
		n_mach(j)=K(j)/dim_mach;
		age0=agemax+1;
		Amax=agemax+1;
		if (n_mach(j) <= 0)
			cerr << "\n ERRORE Inizializzazione!!!" << endl;
		while (n_mach(j) > 0)
		{
			i++;	// Solo le macchine di prima generazione (tt=1)		 
			if (i > N1)
				i=1;
			if (fornit(j) != i)
			{
				g[0][i-1][j-1]++;	// hanno frequenza diversa da zero.	
				gtemp[0][i-1][j-1]++;
				g_c[0][i-1][j-1]++;
				age[0][i-1][j-1]=age0;
				age0--;
				if (age0 < 1)
					age0=agemax+1;
				n_mach(j)--;
				if (flagTCGEN == 1 && flagTC == 2)
					age[0][i-1][j-1]=1;  // Con progresso tecnico det, evitiamo scrapping x et‚àö‚Ä†
				
			}
		}
		
	}
	
	// Inizializzazione D2(t-1) sotto hp steady state uguale x tutte le imprese
	if (flagC==2)
	{
		if (flagTCGEN==1)
			D20=(w(2)/(A1p(i)*a)*(I(2,1)/dim_mach)*N2r*(1-wu-aliq)+wu*w(2)*LS)/(p2(1,1)-((1-wu-aliq)*w(2)/A2(1)));
		if (flagTCGEN==2 && flagTAX==1)
			D20=((w(2)/(A1p(i)*a)+nu*p1(1,1))*I(2,1)/dim_mach*N2r*(1-wu-aliq)+wu*w(2)*LS)/(p2(1,1)-((1-wu-aliq)*w(2)/A2(1)));
		if (flagTCGEN==2 && flagTAX==2)
			D20=((w(2)/(A1p(i)*a)+nu*p1(1,1))*I(2,1)/dim_mach*N2r*(1-wu)+wu*w(2)*LS)/(p2(1,1)-((1-wu)*w(2)/A2(1)));
	}
	else D20=(w(2)/(A1p(i)*a)*(I(2,1)/dim_mach)*N2r)/(p2(1,1)-w(2)/A2(1));
	
	D2=D20/N2r;
	De=D20/N2r;								// Aspettative adattive miopi
	
	
	Ne=theta*D20/N2r;
	N=theta*D20/N2r;
	Nm=(theta*D20/N2r)*p2(1,1);
	S1=(((I(2,1)/dim_mach)*N2r)/N1r)*p1(1,1);	
	S2=D20/N2r*p2(1,1);
	dNtot2=0;
	
	step=N2/N1;							// Inizializzo matrice matching macchine
	cont=0;
	
	if (N2%N1 != 0)
		cerr << "\n\n ERRORE DI INIZIALIZZAZIONE!!!!!" << endl;
	
	Match=0;
	
	for (i=1; i<=N1; i++)
	{
		cont+=step;
		for (j=0; j<step; j++)
		{
			Match(cont-j,i)=1;
			fornit(cont-j)=i;
		}
	}
	
	if (flagTCGEN==2)						// Se progresso tecnico endogeno e inizio turno, devo lanciare innovazione
	{		
		A1top=A0;							// Inizializzo frontiera endogena
		A1ptop=A0;
		A1front=A0;							// Inizializzo frontiera esogena
		A1pfront=A0;
		A1f=A0;
		A1pf=A0;
		RD=nu*S1(2,1);						// Inizializzo R&D.	
		t=0;
		TECHANGEND();						// procedura fuori dal ciclo temporale
		S1=(((I(2,1)/dim_mach)*N2r)/N1r)*p1(1,1);	
		A1top(2)=A1top(1);
		A1ptop(2)=A1ptop(1);
		A1front(2)=A1front(1);
		A1pfront(2)=A1pfront(1);
	}
	else {	  
		RD=0;								// Azzero R&D e domanda lavoro
		Ld1rd=0;
		LD1rdtot=0;
	} 
	
	tao=1;									// Generazione macchinari
}

//-------------------------------------------------------------------------------------------

void INTFILE(void) // Apre file e scrive intestazioni dei file 
{
	ofstream inv_output1(nomefile1);	// Creao file per append successivo
	ofstream inv_output(nomefile);
	ofstream inv_prod1(nomefile5);
	ofstream inv_prod2(nomefile4);
	ofstream inv_prodall1(nomefile6);
	ofstream inv_prodall2(nomefile7);

	inv_output.setf(ios::fixed);
	inv_output.precision(3);
	inv_output.setf(ios::right);
	inv_output.width(6);
	inv_output <<"t";
	inv_output.width(30);
	inv_output <<"Q1+Q2";
	inv_output.width(30);
	inv_output <<"GDP";
	inv_output.width(30);
	inv_output << "EI";
	inv_output.width(30);
	inv_output << "SI";
	//inv_output.width(20);
	//inv_output << "N";
	inv_output.width(25);
	inv_output << "dN";
	inv_output.width(30);
	inv_output << "Q1";
	inv_output.width(30);
	inv_output << "Q2";
	inv_output.width(30);
	inv_output << "C/cpi";
	inv_output.width(30);
	inv_output << "I/ppi";
	// inv_output.width(20);
	// inv_output << "Saldo";
	inv_output.width(35);
	inv_output << "Debtot1";
	inv_output.width(35);
	inv_output << "Debtot2";
	inv_output.width(35);
	inv_output << "Pi1";
	inv_output.width(35);
	inv_output << "Pi2";
	inv_output.width(35);
	inv_output << "NW1";
	inv_output.width(35);
	inv_output << "NW2";
	inv_output.width(20);
	inv_output << "LD (E)";
	inv_output.width(15);
	inv_output << "u";
	inv_output.width(45);
	inv_output << "w";
	inv_output.width(15);
	inv_output << "d_w";
	//inv_output.width(15);
	//inv_output << "ppi";
	inv_output.width(45);
	inv_output << "cpi";
	inv_output.width(15);
	inv_output << "d_cpi";
	inv_output.width(20);
	inv_output << "Am";
	inv_output.width(20);
	inv_output << "A_sd";
	//inv_output.width(20);
	//inv_output << "Am 1";
	//inv_output.width(20);
	//inv_output << "Am 2";
	inv_output.width(15);
	inv_output << "Nexit1";
	inv_output.width(15);
	inv_output << "Nexit2";
	inv_output.width(25);
	inv_output << "w/p";
	inv_output.width(15);
	inv_output << "mi";
	inv_output.width(50);
	inv_output << "G";
	inv_output.width(50);
	inv_output << "Tax";
	inv_output.width(50);
	inv_output << "Def";
	inv_output.width(50);
	inv_output << "Deb";
	inv_output.width(35);
	inv_output << "Def/GDP";
	inv_output.width(35);
	inv_output << "Deb/GDP";
	inv_output.width(15);
	inv_output << "Herf1";
	inv_output.width(15);
	inv_output << "Herf2";
	inv_output.width(15);
	inv_output << "N2lowI";
	inv_output.width(15);
	inv_output << "N2highI";
	inv_output.width(15);
	inv_output << "N2lowEI";
	inv_output.width(15);
	inv_output << "N2ighhEI";
	inv_output.width(15);
	inv_output << "DF1";
	inv_output.width(15);
	inv_output << "DF2" << endl;
	inv_output.close();
	
}

//-------------------------------------------------------------------------------------------

void MACH(void)	// Imprese consumption-good ricevono macchine e aggiornano quantit‚àö‚Ä†. Se non hanno soldi muoiono					

{			
	if (flagTCGEN==1 && flagTC==2 && t==2)	// con progr. tecn. det. bisogna modficare
	{										// modello in modo che al tempo 2 la produttivit‚àö‚Ä† cresca
		cout << "\n flagTC=2 in Mach" << endl;
		for (i=1; i<=N1; i++)
		{
			A(1,i)=A0*(1+uu);
			C(1,i)=w(2)/(A0*(1+uu));
		}
	}
	
	N2eff=N2r;													// Numero imprese effettivo
	for (j=1; j>=N2; j++)										// Ke pu‚àö‚â§ essere > 0 solo se entry nel settore 2
	{															// non ‚àö¬Æ random
		if (Ke(j)>0 && flagENTRY<2)
			cerr << "\n ERRORE entry settore 2 Ke!!!!" << endl;
		if(Ke(j)>0)												// Se impresa ‚àö¬Æ appena entrata non ha K e non ‚àö¬Æ
			N2eff--;											// calcolata nel numero imprese effettivo
	}
	
	if (N2eff==0)
	{
		cerr << "\n ERRORE: nessun sopravvissuto settore 2!!!" << endl;
	}
	
	flag=0;
	p1m=0;
	c1m=0;
	S1m=0;
	
	for (i=1; i<=N1; i++)								// Costi unitari macchine esistenti
	{
		if (flagTCGEN!=2 && RD(2,i) !=0)
			cerr << "\n ERRORE: R&D > 0!!!";		
		S1(1,i)=0;
 		for (tt=t0; tt<=t; tt++)				
		{
			for (j=1; j<=N2; j++)
			{
				if (A(tt,i) > 0)
					C(tt,i)=w(2)/A(tt,i);
				else
				{
					A(tt,i)=A(tt-1,i);
					C(tt,i)=C(tt-1,i);
					cerr << "\n\n A(tt,i) = 0!!!" << endl;
				}
			}	
		}
		if (c1(i) > 0)
			c1(i)=w(2)/(A1p(i)*a);
		else
		{
			A1p(i)=0.01;
			cerr << "\n\n A1p(i) = 0!!!" << endl;
		}
		p1(1,i)=(1+mi1)*c1(i);
		if (p1(1,i) < pmin)
		{
			contppi1++;
			p1(1,i)=pmin;
			if (p1(2,i) < pmin)
				p1(2,i)=pmin;
		}
		p1m+=p1(2,i);
		c1m+=c1(i);
	}

	c1m/=N1;
	p1m/=N1;

	A1m=A1.Sum()/N1r;

	if (flagbug == 1)
	{
		debug << "\n\n\n T E M P O:   " << t << endl;
		debug << "----------------------------------------------------------------" << endl;
		debug << "\n Consegna maccchine..." << endl;
		debug << "\n A:\n\n" << A << endl;
		debug << "\n C:\n\n" << C << endl;
		debug << "\n w(2): " << w(2) << endl;
		debug << "\n c1: " << c1 << endl;
		debug << "\n p1: " << p1 << endl;
		debug << "\n A1p: " << A1p << endl;
		debug << "\n A1: " << A1 << endl;
		for (j=1; j<=N2; j++)								 
		{	
			mat2d=0;
			debug << "\n g   Firm " << j << endl;
			for (i=1; i<=N1; i++)							 
			{	
		 		for (tt=t0; tt<=t; tt++)
				{
					mat2d(tt,i)=g[tt-1][i-1][j-1]; 
				}
			}

			debug << mat2d << endl;
			debug << "\n gtemp   Firm " << j << endl;
			mat2d=0;
			for (i=1; i<=N1; i++)							 
			{												
		 		for (tt=t0; tt<=t; tt++)
				{
					mat2d(tt,i)=gtemp[tt-1][i-1][j-1]; 
				}
			}
			debug << mat2d << endl;
			debug << "\n age   Firm " << j << endl;
			mat2d=0;
			for (i=1; i<=N1; i++)							 
			{												
		 		for (tt=t0; tt<=t; tt++)
				{
					mat2d(tt,i)=age[tt-1][i-1][j-1]; 
				}
			}
			debug << mat2d << endl;
		}
	}

	for (j=1; j<=N2; j++)								// Imprese ricevono macchine, 
	{													// aggiornano stock di capitale e 
		flag=0;
		K(j)+=EI(2,j);									// scrappano macchine che hanno costo produzione superiore
		for (i=1; i<=N1; i++)							// costo produzione superiore 
		{												// prezzo del periodo precedente
		 	for (tt=t0; tt<=t; tt++)
			{
				g[tt-1][i-1][j-1]=gtemp[tt-1][i-1][j-1]; 
				g_c[tt-1][i-1][j-1]=g[tt-1][i-1][j-1];
	
				if (g[tt-1][i-1][j-1] > 0 && w(3)/A(tt,i) > p2(2,j) && flagPC == 1)
				{
					cerr << "\n ERRORE: c(j) < p2(2,j)!!!" << endl;
					K(j)-=g[tt-1][i-1][j-1]*dim_mach;
					g[tt-1][i-1][j-1]=0;
					g_c[tt-1][i-1][j-1]=0;
					gtemp[tt-1][i-1][j-1]=0;
					if (K(j) <= 0 && Ke(j)<=0)
						flag=1;
				}
			}
		}
	} 

	if (flag==1)
		cerr << "\n\n K <= 0!" << endl;

	for (j=1; j<=N2; j++)								
	{
		if (Ke(j) == 0)										// Se impresa ‚àö¬Æ incumbent,
		{													// calcolo nuova produttivit‚àö‚â§, costi
			n_mach(j)=K(j)/dim_mach;						// calcolo nuovo prezzo vendita e
			c2(1,j)=0;										// costo totale macchine acquistate
			A2(j)=0;
			for (i=1; i<=N1; i++)				
			{
				for (tt=t0; tt<=t; tt++)
				{
					if (n_mach(j)>0)
					{
						c2(1,j)+=C(tt,i)*g[tt-1][i-1][j-1]/n_mach(j);
						A2(j)+=A(tt,i)*g[tt-1][i-1][j-1]/n_mach(j);
					}
					else {
						cerr << "\n\n ERRORE nmach = 0!!! " << endl;
					}
				}
			} 
			if (f2(3,j)>0)
				mu2(1,j)=mu2(2,j)*(1+deltami2*((f2(2,j)-f2(3,j))/f2(3,j)));
			else mu2(1,j)=mu2(2,j);
			if (mu2(1,j) <= 0 || mu2(1,j) >=1) {
				cerr << "\n ERRORE: mark-up fuori range!!! " << endl;
				cout << " m2(t): " << mu2(1,j) << "\t\t m2(t-1) " << mu2(2,j) << endl;
				cout << " f2(t-1): " << f2(2,j) << "\t\t f2(t-2): " << f2(3,j) << endl;
			}			
			p2(1,j)=(1+mu2(1,j))*c2(1,j);						// Nuovo prezzo
			if (p2(1,j) < pmin)
			{
				contcpi1++;
				p2(1,j)=pmin;
				if (p2(2,j) < pmin)
					p2(2,j)=pmin;
			}
		}
	}

	if (flagbug == 1)
	{
		debug << "\n EI:\n\n" << EI << endl;
		debug << "\n K:\n\n" << K << endl;
		debug << "\n Ke:\n\n" << Ke << endl;
		debug << "\n n_mach:\n\n" << n_mach << endl;
		debug << "\n c2:\n\n" << c2 << endl;
		debug << "\n A2:\n\n" << A2 << endl;
		debug << "\n p2:\n\n" << p2 << endl;
		debug << "\n mu2:\n\n" << mu2 << endl;
	}

	for (j=1; j<=N2; j++)
	{
		indforn=int(fornit(j));									// Indice fornitore
		if (I(2,j) > 0)
		{
			Cmach(j)=p1(2,indforn)*I(2,j)/dim_mach;				// Costo totale macchine
			S1m+=I(2,j)/dim_mach;
		}
		else Cmach(j)=0;
		Debmax2(j)=phi1*S2(2,j);								// Limite debito imprese
		if (phi1*W2(2,j) > Debmax2(j))
			Debmax2(j)=phi1*W2(2,j);
		if (Debmax2(j) > Deb2(2,j))
			debres2(j)=Debmax2(j) - Deb2(2,j);
		else debres2(j)=0;
		if (W2(2,j) >= Cmach(j))							// Ricchezza >= costo macchine
		{
			W2(1,j)=W2(2,j)-Cmach(j);
			Deb2(1,j)=Deb2(2,j);
			if (indforn > 0)
				S1(1,indforn)+=Cmach(j);
		}
		else												// Ricchezza < costo macchine
		{ 
			if (Cmach(j) <= W2(2,j) + debres2(j))			// Nuovo debito permette pagare 	
			{												// macchine
				Deb2(1,j)=Deb2(2,j)+Cmach(j)-W2(2,j);
				debres2(j)-=Cmach(j)-W2(2,j);
				W2(1,j)=1;
				if (indforn > 0)
					S1(1,indforn)+=Cmach(j);
			}
			else											// Debito non sufficiente a pagare
			{	
				if (indforn > 0)							// macchine
				{
					if (flagdie==0)
						S1(1,indforn)+=W2(2,j)+debres2(j);// Impresa beni capitali subisce perdite
					else
						S1(1,indforn)+=Cmach(j);
				}
				//contQI++;
				W2(1,j)=0;				
				f2(1,j)=0;									// Impresa manifatturiera muore
				f2(2,j)=0;
				f2(3,j)=0;
				debres2(j)=0;
			}
		}
		if (Deb2(1,j) > Debmax2(j) && Deb2(1,j) > Deb2(2,j))
		{
			cerr << "\n ERRORE: Deb > Debmax!!!" << endl;
			cerr << " Deb2(1,j): " << Deb2(1,j) << endl;
			cerr << " Deb2(2,j): " << Deb2(2,j) <<endl;
			cerr << " Debmax2: " << Debmax2(j) << endl;
			cerr << " debres: " << debres2(j) << endl;
			cerr << " W2(1,j): " << W2(1,j) << endl;
			cerr << " W2(2,j): " << W2(2,j) << endl;
		}
	}

	if (flagdie==1)
	{
		for (i=1; i<=N1; i++)
			S1(1,i)+=die(i);  // Ricevono pagamento imprese morte fine periodo precedente
	}

	if (flagbug == 1)
	{
		debug << "\n Pagamento macchine..." << endl;
		debug << "\n Cmach:\n\n" << Cmach << endl;
		debug << "\n Debmax:\n\n" << Debmax2 << endl;
		debug << "\n S1:\n\n" << S1 << endl;
		debug << "\n RD:\n\n" << RD << endl;
		debug << "\n W2:\n\n" << W2 << endl;
		debug << "\n Deb:\n\n" << Deb2 << endl;
		debug << "\n S2:\n\n" << S2 << endl;
		debug << "\n f2:\n\n" << f2 << endl;
	}

	for (j=1; j<=N2; j++)				
	{
		dt(j)=Deb2(1,j);
		machtemp=n_mach(j);
		contmach=0;
	 	for (i=1; i<=N1; i++)   
		{ 
			for (tt=t0; tt<=t; tt++)
			{
				machtemp-=g[tt-1][i-1][j-1];
				contmach+=g[tt-1][i-1][j-1];
			}
		}
		
		if (machtemp != 0)
		{
			cerr << "\n\n ERRORE: K != g[]!!!   UNO   " << machtemp << endl;
			cout << "\n K: " << K(j) << endl;
			cout << " n_mach: " << n_mach(j) << endl;
			cout << " I: " << I(2,j)/dim_mach << endl;
			cout << "\n Impresa " << j << endl;
			cout << " Fornitore: " << fornit(j) << endl;
			
			//n_mach(j)=contmach;
			//K(j)=n_mach(j)*dim_mach;
			
		}
	}
	for (j=1; j<=N2; j++)
	{
		if (dt(j) != Deb2(1,j))
			cerr << "\n ERRORE: Debito cambiato MACH!!!" << endl;
	}	
}

//------------------------------------------------------------------------------------------

void BROCHURE(void)				// Calcolo competitivit‚àö‚Ä† industria macchinari + invio brochure
{
	ftot=0;	
	flag=0;
	nclient=0;

	for (j=1; j<=N2; j++)							// Se impresa che produce beni di consumo non ha 
	{												// il fornitore, perch‚àö¬Æ ‚àö¬Æ morto nel turno precedente
													// ne ottiene uno in maniera random
		if (fornit(j) < 1 || fornit(j) > N1)
		{
			fornit(j)=int(ran1(p_seed)*10000)%N1+1;			
			Match(j,fornit(j))=1;
		}
	}

	if (flagbug == 1)
	{
		debug <<"\n Brochure: " << endl;
		debug << "\n Fornitore:\n\n" << fornit << endl;
		debug << "\n Match:\n\n" << Match << endl;
	}
	
	for (i=1; i<=N1; i++)		// Per ogni impresa che produce macchinari, determinano quali
	{							// potenziali clienti riceveranno la brochure
		for (j=1; j<=N2; j++)
			nclient(i)+=Match(j,i);								// Numero clienti
		if (nclient(i) <=0)
			flag=1;
		//newbroch=int(nclient(i) + ROUND(nclient(i)*Gamma));		// Numero potenziali clienti
		newbroch=int(ROUND(nclient(i)*Gamma));
		if (newbroch==0)
			newbroch++;
		//if (nclient(i) + newbroch >= N2)						// Limite max brochure
		//	newbroch=N2-nclient(i);
		if (f1(2,i)>f1max)										// Se f1max < 1 c'è antitrust
			newbroch=0;
		while (newbroch > 0)
		{
			rni=int(ran1(p_seed)*10000)%N2+1;
			Match(rni,i)=1;
			newbroch--;
		} 
	} 
	if (flag==1)
		cerr << "\n\n ERRORE: Numero clienti pari a zero!!!" << endl;

	if (flagbug == 1)
	{
		debug << "\n Dopo invio brochure..." << endl;
		debug << "\n nclient:\n\n" << nclient << endl;
		debug << "\n Match:\n\n" << Match << endl;
		debug << "\n A1:\n\n" << A1 << endl;
	}

	for (j=1; j<=N2; j++)				// Ogni impresa che produce beni di consumo sceglie
	{									// fornitore di fiducia
		indforn=int(fornit(j));
		for (i=1; i<=N1; i++)
		{
			if (A1(i) > 0)
			{
			//if (Match(j,i)==1 && A1(i)/p1(1,i) > A1(indforn)/p1(1,indforn)) // Impresa sceglie fornitore con max A/p
				if (Match(j,i)==1 && p1(1,i)*(w(2)/A1(i))*b < p1(1,indforn)*(w(2)/A1(indforn))*b) //Impresa sceglie fornitore con minore p*c*b 
					indforn=i;
			}
			else cerr << "\n\n ERRORE: A1(i) = 0!!!" << endl;
		}
		fornit(j)=indforn;
		for (i=1; i<=N1; i++)
		{
			if (i != indforn)
				Match(j,i)=0;
		}
	}

	nclient=0;
	for (i=1; i<=N1; i++)
	{
		for (j=1; j<=N2; j++)
		{
			if (Match(j,i) == 1)
				nclient(i)++;
		}
	}

	if (flagbug == 1)
	{
		debug << "\n Dopo matching, clienti cambiano..." << endl;
		debug << "\n Fornitore:\n\n" << fornit << endl;
		debug << "\n Match:\n\n" << Match << endl;
		debug << "\n nclient:\n\n" << nclient << endl;
	}

	for (j=1; j<=N2; j++)
	{
		if (dt(j) != Deb2(1,j))
			cerr << "\n ERRORE: Debito cambiato BROCH!!!" << endl;
	}
	
	nclientmax=0;
	for (i=1; i<=N1; i++) // calcolo numero di clienti massimi che ha impresa beni capitali
	{
		if (nclient(i) > nclientmax)
			nclientmax=nclient(i);		
	}
	if (nclientmax/N2r >= 0.99)
	{
		//cout << "\n Quota mercato impresa leader beni capitali: " << nclientmax/N2r;
		contMONOP++;
	}
}

//------------------------------------------------------------------------------------------

void INVEST(void)									// Investimenti

{ 
	//cout <<"\nInvest..." << endl;
	EId=0;
	flag=0;
	agetop=0;
	
	for (j=1; j<=N2; j++)
	{
		if (dt(j) != Deb2(1,j))
			cerr << "\n ERRORE: Debito cambiato INV1!!!" << endl;
	}

	npayback=0;

	for (j=1; j<=N2; j++)
	{
		if (dt(j) != Deb2(1,j))
		{
			if (echo==1)
				cerr << "\n BUG INV!!!" << endl;
			//cout << " dt: " << dt(j) << endl;
			//cout << " Deb2(1,j): " << Deb2(1,j) << endl;
			//cout << " Deb2(2,j): " << Deb2(2,j) << endl;
			Deb2(1,j)=dt(j);
		}

		EXPECT();
		if (Kexp(j)==1 && flagENTRY >=2){			// Se impresa entra periodo t-2
			De(1,j)=u*K(j);							// De = u*K
			if (De(1,j)<=0)
				cerr << "\n\n ERRORE: Kexp=1 & De=0!!! " << De(1,j) << K(j) << j << endl;
		}
		Ne(j)=De(1,j)*theta;						// scorte attese
		
		Qd(j)=De(1,j)+Ne(j)-N(2,j);					// produzione desiderata
		
		if (Qd(j) < 0) 
			Qd(j)=0;

		if (t==1)
			Kd(j)=Qd(j);							// capitale desiderato
		else
			Kd(j)=Qd(j)/u;
		
		Ktrig(j)=K(j)*(1+alfa);						// calcolo soglia per investimento
		Ktrig(j)=ROUND(Ktrig(j)/dim_mach)*dim_mach;	// arrotondamento K trigger
		Ktop(j)=K(j)*(1+alfasup);					// calcolo limite superiore per investimento
		Ktop(j)=ROUND(Ktop(j)/dim_mach)*dim_mach;	// arrotondamento K top
				
		if (Kd(j) >= Ktrig(j))						// calcolo investimento espansione 
		{
			if (alfa > 0)
				EId(j)=Ktrig(j)-K(j);				// desiderato
			else									// se non c'‚àö¬Æ trigger, comunque gli 
			{										// investimenti incontrano un limiter superiore
				if (alfasup > 0 && Kd(j) > Ktop(j))
				{
					EId(j)=Ktop(j)-K(j);
					contalfasup++;
				}
				else
					EId(j)=floor((Kd(j)-K(j))/dim_mach)*dim_mach;
			}
		}
		else EId(j)=0;
		
		if (flagENTRY>=2 && Ke(j)>0)				// Se impresa ‚àö¬Æ entrata periodo t-1
		{											// si deve comprare il K
			Kd(j)=Ke(j);
			EId(j)=Kd(j);
		}
		
		if (Qd(j) > K(j))							// Se produzione desiderata > capacit‚àö‚Ä†
		{											// produttiva dopo investimento,
			Qd(j)=K(j);								// va ridotta!
		}

		SCRAPPING();								// calcolo investimenti sostituzione

		if (Qd(j) > 0 && Qd(j) < K(j))				// Produzione minore capitale: uso
		{
			COSTPROD();								// macchine pi‚àöœÄ produttive
		}
		else
		{
			A2e(j)=A2(j);
			c2e(j)=c2(1,j);

		}

		PRODORD();						// calcolo investimenti totali domandati con credito
	}

	if (flagbug == 1)
	{
		debug << "\n Produzione & Investimenti" << endl;
		debug << "\n Kexp:\n\n" << Kexp << endl;
		debug << "\n Ke:\n\n" << Ke << endl;
		debug << "\n K:\n\n" << K << endl;
		debug << "\n Kd:\n\n" << Kd << endl;
		debug << "\n Qd:\n\n" << Qd << endl;
		debug << "\n EId:\n\n" << EId << endl;
		debug << "\n SId:\n\n" << SId << endl;
		debug << "\n A2e:\n\n" << A2e << endl;
		debug << "\n c2e:\n\n" << c2e << endl;
		debug << "\n A2:\n\n" << A2 << endl;
		debug << "\n c2:\n\n" << c2 << endl;
		debug << "\n Q2:\n\n" << Q2 << endl;
		debug << "\n EI:\n\n" << EI << endl;
		debug << "\n SI:\n\n" << SI << endl;
		debug << "\n I:\n\n" << I << endl;
		debug << "\n W2:\n\n" << W2 << endl;
		debug << "\n Deb:\n\n" << Deb2 << endl;
	}
	if (flag==1 && echo==1)
		cout << "\n I < EId + SId!!!" << endl;
	if (flag==2 && echo==1)
		cout << "\n Q2 < Qd!!!" << endl;

	Itot=0;
	ninv=0;
	for (j=1; j<=N2; j++)
	{
		Itot+=I(1,j);
		if (I(1,j) < 0 || EI(1,j) < 0 || SI(1,j) < 0)
		{
			cerr <<"\n ERRORE: I < 0 " << endl;
		}
		if (I(1,j)==0)
			ninv++;
	}
	Itot/=dim_mach;

	if (agetop > agemax+1 && echo==1)
		cerr << "\n Age scrapping macchinari: " << agetop << endl;
}

//------------------------------------------------------------------------------------------

void EXPECT(void)						// Aspettative domanda imprese industria 2

{
	switch( flagEXP ) 
	{
    case 0:							// Aspettative miopi
        De(1,j) = D2(2,j);
        break;
    case 1 :						// Aspettative miopi con memoria 1-4
        De(1,j) = beta1*D2(2,j) + beta2*D2(3,j) + beta3*D2(4,j) + beta3*D2(5,j);
        break;
    case 2 :						// Aspettative accelerative GD 
		if (D2(3,j)<=0)
			D2(3,j)=1;
		De(1,j) = (1 + beta*(D2(2,j)-D2(3,j))/D2(3,j)) * D2(2,j);
        break;
	case 3 :						// Aspettative adattive
		De(1,j) = De(2,j) + delta1*(D2(2,j) - De(2,j));
        break;
	case 4 :						// Aspettative estrapolative - accelerative
		if (D2(3,j)<=0)
			D2(3,j)=1;
		if (GDP(3)<=0)
			GDP(3)=1;
		De(1,j) = (1 + delta5*(D2(2,j) - D2(3,j))/D2(3,j) + delta6*(GDP(2) - GDP(3))/GDP(3))*D2(2,j);
		break;
	default : 
		cerr << "\n\n ERRORE aspettative!" << endl;
	}
}

//------------------------------------------------------------------------------------------

void SCRAPPING(void)							// Calcola investimenti di rimpiazzo desiderati

{	 
 	SId(j)=0;
	indforn=int(fornit(j));
	
	for (i=1; i<=N1; i++)
	{	 
		for (tt=t0; tt<=t; tt++)
		{ 
			if (age[tt-1][i-1][j-1] > agetop)
				agetop=age[tt-1][i-1][j-1];
			C_pb[tt-1][i-1][j-1]=0;
			g_pb[tt-1][i-1][j-1]=0;								// Se macchina non deve essere
			if (g[tt-1][i-1][j-1] > 0 && A(tt,i) < A1(indforn))	// sostituita (>b) g_pb rimane 0
			{													// altrimenti memorizzo numero 
				//payback=p1/(w(2)/A(tt,i)-w(2)/A1(indforn));	// macchine da sostituire. 
				if (w(2) > 0 && A(tt,i) > 0 && A1(indforn) > 0)
					payback=p1(1,indforn)/(w(2)/A(tt,i)-w(2)/A1(indforn));
				else cerr << "\n\n ERRORE: payback divisione per zero!!!" << endl;
				
				if (payback <= b)							 							
				{												
					npayback++;
					if (contSI==0 && echo==1)
						cout << "\n\n Primo Scrapping!!!" << endl;
					contSI++;
					g_pb[tt-1][i-1][j-1]=g[tt-1][i-1][j-1];	
		 			C_pb[tt-1][i-1][j-1]=C(tt,i);			// costo produzione macchine sostituire
					SId(j)+=dim_mach*g_pb[tt-1][i-1][j-1];	// Invest. sostituz. desiderati					
				}
			}
			if (g[tt-1][i-1][j-1] > 0 && g_pb[tt-1][i-1][j-1]==0 && age[tt-1][i-1][j-1] > agemax)
			{ // Macchine scrappate anche se superano et‚àö‚Ä† massima
					//npayback++;
				g_pb[tt-1][i-1][j-1]=g[tt-1][i-1][j-1];	
	 			C_pb[tt-1][i-1][j-1]=C(tt,i);			// costo produzione macchine sostituire
				SId(j)+=dim_mach*g_pb[tt-1][i-1][j-1];	// Invest. sostituz. desiderati
			}
		}
	}
}

//-------------------------------------------------------------------------------------------

void COSTPROD(void)  // Calcola costo unitario investimenti di rimpiazzo desiderati
{
	A2e(j)=0;
	c2e(j)=0;
	nmachprod=ceil(Qd(j)/dim_mach);			// Numero di macchine necessarie per produzione
	nmp_temp=nmachprod;
	
	machtemp=K(j)/dim_mach;
	contmach=0;
	for (i=1; i<=N1; i++)   
	{ 
		for (tt=t0; tt<=t; tt++)
		{
			machtemp-=g_c[tt-1][i-1][j-1];
			contmach+=g_c[tt-1][i-1][j-1];
		}
	}
	if (machtemp != 0 || K(j)==0)
	{
		cerr << "\n\n ERRORE: K != g[]!!!   " << machtemp << endl; 
		//n_mach(j)=contmach;
		//K(j)=n_mach(j)*dim_mach;
	}

	while (nmp_temp > 0)				// Cicla finch‚àö¬Æ non troviamo il mix di macchine
	{									// pi‚àöœÄ produttivo
		Amax=0;
		imax=0;
		jmax=0;
		tmax=0;
		
		for (i=1; i<=N1; i++)			// Trovo macchina con produttivit‚àö‚Ä† massima
		{
			for (tt=t0; tt<=t; tt++)
			{
				if (g_c[tt-1][i-1][j-1] > 0 && A(tt,i) > Amax)
				{							
					Amax=A(tt,i);
					imax=i;
					jmax=j;
					tmax=tt;
				}
			}
		}

		if (nmachprod>0)
		{
			if (g_c[tmax-1][imax-1][jmax-1] >= nmp_temp)	// Produzione resid. <= dim. macchina
			{
				A2e(j)+=A(tmax,imax)*nmp_temp/nmachprod;
				c2e(j)+=C(tmax,imax)*nmp_temp/nmachprod;
				g_c[tmax-1][imax-1][jmax-1]-= nmp_temp;
				nmp_temp=0;									// usciamo da while
			}
			else											// Produzione resid. > dim. macchina
			{
				A2e(j)+=A(tmax,imax)*g_c[tmax-1][imax-1][jmax-1]/nmachprod;
				c2e(j)+=C(tmax,imax)*g_c[tmax-1][imax-1][jmax-1]/nmachprod;
				nmp_temp-=g_c[tmax-1][imax-1][jmax-1];
				g_c[tmax-1][imax-1][jmax-1]=0;				// Eliminiamo macchina + produttiva
			}												// in modo che non sia scelta
		}
		else cerr << "\n\n ERRORE: nmachprod = 0!!!" << endl;
	}
}

//-------------------------------------------------------------------------------------------

void PRODORD(void)		// Calcola produzione effettiva e ordini investimenti				
{	
 	EI(1,j)=0;		// Nota: l'impresa ha uscite effettive per la produzione, mentre 					
 	SI(1,j)=0;		// per gli investimenti effettua calcoli prudenziali				
	I(1,j)=0;

	if (Deb2(1,j) > Debmax2(j) && Deb2(1,j) > Deb2(2,j))
	{
		cerr << "\n ERRORE: Deb > Debmax!!! +++" << endl;
	}

	if (c2e(j)*Qd(j) <= W2(1,j))	// Determina produzione
	{
		Q2(1,j)=Qd(j);
		W2(1,j)-=c2e(j)*Qd(j);
	}
	else
	{
		if (c2e(j)*Qd(j) <= W2(1,j) + debres2(j))		 	
		{											
			Deb2(1,j)+=c2e(j)*Qd(j)-W2(1,j);
			debres2(j)-=c2e(j)*Qd(j)-W2(1,j);
			Q2(1,j)=Qd(j);
			W2(1,j)=1;
		}
		else											
		{
			if (c2e(j) > 0)
				Q2(1,j)=(W2(1,j) + debres2(j))/c2e(j);
			else cerr << "\n\n c2e(j) == 0!!!" << endl;
			W2(1,j)=1;	
			Deb2(1,i)+=debres2(j);
			debres2(j)=0;
			if (Q2(1,j) < 1)		// Se impresa non riesce a produrre, muore
			{
				f2(1,j)=0;									
				f2(2,j)=0;
				f2(3,j)=0;
				W2(1,j)=0;	
			}
			flag=2;
			if (I(2,j) > 0)  // conto quante imprese che hanno pagate macchine non riescono a produrre
				contQI++;
			if (flagENTRY>=2 && Kexp(j)==1)
			{
				cerr << "\n\n ERRORE: Kexp=1 e Debito x Prod.!!!" << endl;
			}
		}
	}

	if (flagENTRY>=2 && Ke(j)>0)
		Ld2(j)=0;
	else
	{
		if (A2e(j) > 0)
			Ld2(j)=Q2(1,j)/A2e(j);				// Domanda lavoro
		else cerr << "\n\n A2e(j) > 0" << endl;
	}
	
	NW=W2(1,j);								// Calcoli teorici e prudenziali per determinare investimenti
	prestmax=phi2*mol(j);					// Limite max a nuovi prestiti (ratio: prudenza)
	if (prestmax < 0)
		prestmax=0;
	if (debres2(j) > prestmax)
		prestmax=debres2(j);
	indforn=int(fornit(j));
	if ((EId(j)/dim_mach)*p1(1,indforn) < NW)
	{
		EI(1,j)=EId(j);
		NW-=(EId(j)/dim_mach)*p1(1,indforn);
	}
	else 
	{	
		p1prova=p1(1,indforn);				// Necessario altrimenti mi crea function p1!
		if (p1prova <=0)
			cerr << "\n\n ERRORE p1 <=0!!!" << endl;
		if ((EId(j)/dim_mach)*p1prova-NW <= prestmax)
		{
			EI(1,j)=EId(j);
			prestmax-=(EId(j)/dim_mach)*p1(1,indforn)-NW;
			NW=0;
	    }
	    else 
		{
			EI(1,j)=floor((NW + prestmax)/p1(1,indforn))*dim_mach;
			if (flagENTRY>=2 && Ke(j)>0 && EI(1,j)==0)		// Impresa appena entrata..
			{
				if (NW+prestmax>=p1(1,indforn))				// ... compra 1 macchina
					EI(1,j)=dim_mach;
				else										// ... o muore
				{					
					Ke(j)=0;
					f2(1,j)=0;
					f2(2,j)=0;
					f2(3,j)=0;
					cerr << "\n ERRORE: imprese entrate non riescono ad investire!!!" << endl;
				}
			}
			prestmax=0;
			NW=0;
		}
	}
	
	if (flagENTRY>=2 && Ke(j)>0 && EI(1,j)<=0) {
		cerr << "\n ERRORE: imprese entrate non investono!!!" << endl;
	}
	
	indforn=int(fornit(j));
	p1prova=p1(1,indforn);					// Necessario altrimenti mi crea function p1!
	if (p1prova <=0)
		cerr << "\n\n ERRORE p1 <=0!!!" << endl;
	if ( (SId(j)/dim_mach)*p1prova < NW)
	{
		SI(1,j)=SId(j);
	}
	else 
	{			
		if ((SId(j)/dim_mach)*p1prova-NW <= prestmax)
		{
			SI(1,j)=SId(j);
		}
		else 
		{
			SI(1,j)=floor((NW + prestmax)/p1(1,indforn))*dim_mach;
		}
	}

	I(1,j)=EI(1,j)+SI(1,j);			// Investimenti domandati totali impresa
	if (flag!=2)
	{
		if (EI(1,j) < EId(j) || SI(1,j) < SId(j))
			flag=1;
	}
}

//----------------------------------------------------------------------------------------

void PRODMACH(void)						// Produzione imprese macchinari
{
	Q1=0;
	D1=0;
	flag=0;
	Itot=0;
	D1tot=0;

	for (i=1; i<=N1; i++)				
	{
		if (W1(1,i) <= 0)
			cerr << "\n\n ERRORE: W1 < 0 prima Prodmach " << W1(1,i);

		Debmax1(i)=phi1*S1(2,i);
		if (phi1*W1(2,i) > Debmax1(i))
			Debmax1(i)=phi1*W1(2,i);
		if (Debmax1(i) > Deb1(2,i))
			debres1(i)=Debmax1(i) - Deb1(2,i);
		else debres1(i)=0;
		for (j=1; j<=N2; j++)					// Calcolo domanda di ogni impresa
		{
			Itot+=I(1,j);
			if (Match(j,i) == 1)
			{
				D1(i)+=I(1,j)/dim_mach;
				D1tot+=I(1,j)/dim_mach;
			}
		}
		if (c1(i)*D1(i) <= W1(2,i)+S1(1,i)-RD(2,i))		// Data situazione finanziaria, calcolo
		{												// produzione imprese macchinari
			Q1(i)=D1(i);
			W1(1,i)=W1(2,i)-c1(i)*D1(i)+S1(1,i)-RD(2,i);
			Deb1(1,i)=Deb1(2,i);
			
		}
		else
		{
			if (c1(i)*D1(i) <= W1(2,i)+S1(1,i)-RD(2,i)+debres1(i))
			{
				Q1(i)=D1(i);
				Deb1(1,i)+=c1(i)*D1(i)-W1(2,i)-S1(1,i)+RD(2,i);
				debres1(i)-=c1(i)*D1(i)-W1(2,i)-S1(1,i)+RD(2,i);
				W1(1,i)=1;

			}
			else 
			{
				if (flagdieP==0)
				{
					if (c1(i)>0)
						Q1(i)=floor((W1(2,i)+S1(1,i)-RD(2,i)+debres1(i))/c1(i));	// Credit constrained
					else cerr << "\n\n c1(i) = 0 !!! " << endl;
				}
				else
					Q1(i)=D1(i);										// No credit constrain
				Deb1(1,i)+=c1(i)*D1(i)-W1(2,i)-S1(1,i)+RD(2,i);		
				W1(1,i)=1;

				if (Q1(i) >= 1)
				{
					//Deb1(1,i)+=debres1(i);
					debres1(i)=0;
					W1(1,i)=1;
				}
				else								// Impresa muore
				{
					Q1(i)=0;
					if (D1(i)>0)
					{
						W1(1,i)=0;
						W1(2,i)=0;
					}
				}
			}
		}

		if (D1(i) > Q1(i))							// Domanda > produzione
		{
			flag=1;
			for (j=1; j<=N2; j++)
			{										// Imprese beni di consumo sono razionate
				if (Match(j,i) == 1)
				{
					if (D1(i) > 0)
						I(1,j)=floor((I(1,j)/dim_mach)*Q1(i)/D1(i))*dim_mach;
					else cerr << "\n\n D1(i) = 0!!!" << endl;
					if (I(1,j) < EI(1,j))			// Dato razionamento, ripartizione investimenti
						EI(1,j)=I(1,j);
					SI(1,j)=I(1,j)-EI(1,j);
					if (flagENTRY>=2 && Ke(j)>0 && EI(1,j)==0)	// se impresa appena entrata si compra
						EI(1,j)=dim_mach;						// 1 macchina
				}
			}
		}
		if (A1p(i)>0)
			Ld1(i)=Q1(i)/(A1p(i)*a);					// Domanda lavoro
		else cerr << "\n\n A1p(i)=0!!!" << endl;
		if (W1(1,i) <= 0)
			cerr << "\n\n ERRORE: W1 < 0 dopo Prodmach " << W1(1,i);
	}

	if (flag==1 && echo==1)
		cout << "\n Produzione macchinari razionata!" << endl;

	if (flagbug == 1)
	{
		debug << "\n Produzione macchine..." << endl;
		debug << "\n D1:\n\n" << D1 << endl;
		debug << "\n Q1:\n\n" << Q1 << endl;
		debug << "\n EI:\n\n" << EI << endl;
		debug << "\n SI:\n\n" << SI << endl;
		debug << "\n I:\n\n" << I << endl;
	}

	
	LABOR();								// Calcolo domanda di lavoro nei due settori
	
	if (flagbug == 1)
	{
		debug << "\n Ld2:\n\n " << Ld2 << endl; 
		debug << "\n LD1: " << LD1tot << endl;
		debug << "\n LD2: " << LD2tot << endl;
		debug << "\n LD1rd: " << LD1rdtot;
		debug << "\n LS: " << LS << endl;
	}
	
	Itot/=dim_mach;
	
	for (j=1; j<=N2; j++)				// Aggiorno matrice temporanea frequenze macchine 
	{	
		if (SI(1,j) > 0)				// Cancellazione macchine scrappate
		{
			CANCMACH();
		}
		
		for (i=1; i<=N1; i++)
		{
			for (tt=t0; tt<=t; tt++)
			{
				if (gtemp[tt-1][i-1][j-1] == 0)
					age[tt-1][i-1][j-1]=0;
			}
		}
		indforn=int(fornit(j)); 
		gtemp[t-1][indforn-1][j-1]+=I(1,j)/dim_mach;
		
		if (I(1,j) > 0)
		{
			if (gtemp[t-1][indforn-1][j-1] == I(1,j)/dim_mach)
				age[t-1][indforn-1][j-1]=0;
		}
	}
	
	if (flagbug == 1)
	{
		for (j=1; j<=N2; j++)
		{
			debug << "\n age Firm  " << j << endl;
			for (tt=t0; tt<=t; tt++)				
			{ 
				for (i=1; i<=N1; i++)   
				{ 
					mat2d(tt,i)=age[tt-1][i-1][j-1]; 
				}
			}
			debug << mat2d << endl;
		}
	}
}

//-----------------------------------------------------------------------------------------

void LABOR(void)									// Calcola domanda di lavoro nei 2 settori
{
	LD2tot=Ld2.Sum();
	LD1tot=Ld1.Sum();
	
	if (flagEXTRA==1 && LD1rdtot+LD2tot+LD1tot>LS)
	{
		LS*=1.02;
	}
	LSe=LS;

	if (LD1rdtot < LSe)
		LSe-=LD1rdtot;
	else 
	{
		cerr << "\n\n ERRORE: LD1rdtot > LS!!!   " << LD1rdtot;
	}
	
	if (LD2tot + LD1tot <= LSe)
		LSe=LSe-LD1tot-LD2tot;							// Disoccupazione
	else
	{
		if (echo==1)
			cout << "\n\n Piena occupazione!!!" << endl;
		contFE++;
		
		for (j=1; j<=N2; j++)
		{
			Q2(1,j)=Q2(1,j)*LSe/(LD1tot+LD2tot);						// Riduco produzione imprese
			Ld2(j)=Ld2(j)*LSe/(LD1tot+LD2tot);							// Ricalcolo domanda lavoro	
		}
		
		for (i=1; i<=N1; i++)
		{
			Qpast=Q1(i);
			if (Qpast > 0)
			{
				Q1(i)=floor(Q1(i)*LSe/(LD1tot+LD2tot));
				Ld1(i)=Ld1(i)*LSe/(LD1tot+LD2tot);					// Domanda lavoro
				for (j=1; j<=N2; j++)
				{					
					if (Match(j,i) == 1)
					{
						I(1,j)=floor((I(1,j)/dim_mach)*Q1(i)/Qpast)*dim_mach;
						if (flagENTRY>=2 && Ke(j)>0 && I(1,j)==0)	// E' necessario per evitare che chi entra rimanga senza K
						{
							I(1,j)=dim_mach;
							Q1(i)++;
							Ld1(i)+=1/(A1p(i)*a);					// aggiungo lavoratori per produrre macchina
							//cout << "\n Nuova impresa senza K per arrotondamento! Forzato K=1!" << endl;
						}
						if (I(1,j)<EI(1,j))			// Dato razionamento, ripartizione investimenti
						{
							EI(1,j)=I(1,j);
						}
						SI(1,j)=I(1,j)-EI(1,j);
					}
				}
			}
			//LD2tot=Ld2.Sum();
			//LD1tot=Ld1.Sum();
		}
		LD2tot=Ld2.Sum();
		LD1tot=Ld1.Sum();
	}
	LD=LD1tot+LD2tot+LD1rdtot;				
	LD2=LD1tot+LD2tot;
	if (LD >= LS*(1+0.02))//|| LD < LS-0.1)
	{
		cerr << "\n\n ERRORE: LD > LS dopo rescaling!!!" << "\t\tLD: " << LD << "\t\t LS: " << LS << endl;
	}
	
}



//------------------------------------------------------------------------------------------

void CANCMACH(void)				

 // Elimina dalla matrice delle frequenze assolute macchine effettivamente scrappate.
 // Se l'impresa non ha abbastanza risorse finanziarie per sostituire tutte le 
 // macchine; ci sar‚àö‚Ä† una sostituzione parziale: le macchine saranno sostituite in
 // base al loro costo di produzione.
 
{
 scrapmax=SI(1,j)/dim_mach;			// numero max di macchine che si possono scrappare

 for (i=1; i<=N1; i++)				// Con le due for trovo macchina con costo massimo,
 {									// cio‚àö¬Æ quella che dovr‚àö‚Ä† essere sostituita per prima
	 for (tt=t0; tt<=t; tt++)	
	 {
		 if (scrapmax > 0 && g_pb[tt-1][i-1][j-1] > 0 && age[tt-1][i-1][j-1]>agemax)
		 {
			 if (g_pb[tt-1][i-1][j-1] >= scrapmax)	// Frequenze assoluta macchina con costo 
			 {												// max >= numero macchine che si possono sostituire
 		 		 gtemp[tt-1][i-1][j-1]-=scrapmax;	// Riduco frequenze assolute macchine sostituite
		 		 g_pb[tt-1][i-1][j-1]-=scrapmax;
		 		 scrapmax=0;								// Esco da ciclo while
			 }
	 		 else											// Frequenze assoluta macchina con costo max <
			 {												// numero macchine che si possono sostituire
		 		 scrapmax-=g_pb[tt-1][i-1][j-1];	// Riduco contatore while
 		 		 gtemp[tt-1][i-1][j-1]-=g_pb[tt-1][i-1][j-1];				// Azzero frequenza assoluta macchinario
		 		 g_pb[tt-1][i-1][j-1]=0;
			 }

		 }			
	 }
 }
										// Posizione macchine nella matrice conta
 while (scrapmax > 0)					// sono scrappate prima macchine con indice basso
 {
	 cmax=0;
	 for (i=1; i<=N1; i++)				// Con le due for trovo macchina con costo massimo,
	 {									// cio‚àö√£ quella che dovr‚Äö√Ñ¬∞ essere sostituita per prima
	 	 for (tt=t0; tt<=t; tt++)	
		 {
			 if (g_pb[tt-1][i-1][j-1] > 0 && C_pb[tt-1][i-1][j-1] > cmax)
			 {
				 ind_i=i;					// Memorizzo indici macchina con costo massimo
				 ind_tt=tt;
				 cmax=C_pb[tt-1][i-1][j-1]; // Aggiorno costo massimo
			 }			
		 }
	 }
 	 
	 if (g_pb[ind_tt-1][ind_i-1][j-1] >= scrapmax)	// Frequenze assoluta macchina con costo 
	 {												// max >= numero macchine che si possono sostituire
 		 gtemp[ind_tt-1][ind_i-1][j-1]-=scrapmax;	// Riduco frequenze assolute macchine sostituite
		 g_pb[ind_tt-1][ind_i-1][j-1]-=scrapmax;
		 scrapmax=0;								// Esco da ciclo while
	 }
	 else											// Frequenze assoluta macchina con costo max <
	 {												// numero macchine che si possono sostituire
		 scrapmax-=g_pb[ind_tt-1][ind_i-1][j-1];	// Riduco contatore while
 		 gtemp[ind_tt-1][ind_i-1][j-1]=0;			// Azzero frequenza assoluta macchinario
		 g_pb[ind_tt-1][ind_i-1][j-1]=0;
	 }
 }
 }	

//---------------------------------------------------------------------------------------------

void COMPET2(void)				// Calcolo competitivit‚àö‚Ä† e market share industria manifatturiera 
										
{
	//cout << "\nCompet2..." <<endl;
	
	Em2(1)=0;
	
	if (flagENTRY < 2)								// Random Copy
	{
		p2m=p2.Sum()/N2r;
		l2m=l2.Sum()/N2r;
	}
	else											// Chi entra compra beni capitali
	{
		p2m=0;
		l2m=0;
		N2eff=0;
		for (j=1; j<=N2; j++)
		{
			if (Ke(j) == 0)							// Calcolo media escludendo imprese 
			{										// appena entrate
				p2m+=p2(1,j);
				l2m+=l2(j);
				N2eff++;
			}
		}
		if (N2eff>0)
		{
			p2m/=N2eff;
			l2m/=N2eff;
		}
		else cerr << "\n\n N2eff = 0!!! " << endl;
	}

	ftot=0;
	
	for (j=1; j<=N2; j++)	
	{
		E2(1,j)=-omega1*p2(1,j)/p2m-omega2*l2(j)/l2m;
		ftot(1)+=f2(1,j);
		ftot(2)+=f2(2,j);
		ftot(3)+=f2(3,j);
	}
	
	if (ftot(1)>0)
	{
		for (j=1; j<=N2; j++)			// Riscalo competitivit‚àö‚Ä† tenendo conto imprese morte
		{
			f2(1,j)/=ftot(1);
			f2(2,j)/=ftot(2);
			f2(3,j)/=ftot(3);
			Em2(1)+=E2(1,j)*f2(2,j);
		}
	}
	else 
	{
		cerr << "\n ERRORE: f2 = 0 per tutte le imprese!!!" << endl;
		Em2(1)=Em2(2);
	}
	
	ftot=0;
	if (Em2(1)==0)
		cerr << "\n\n Em2(1)=0!!! " << endl;

	for (j=1; j<=N2; j++)
	{
		if (flagENTRY>=2 && Ke(j)>0)							// Imprese appena entrate non producono
		{														// e hanno market share nulla
			f2(1,j)=0;
			f2(2,j)=0;
			f2(3,j)=0;
		}
		else
		{
			f2(1,j)=f2(2,j)*(1+chi*(E2(1,j)-Em2(1))/Em2(1));	// aggiorno market share
			if (f2(1,j) < exit2)
			{	
				f2(1,j)=0;
				f2(2,j)=0;
				f2(3,j)=0;
			}
		}
		ftot(1)+=f2(1,j);
		ftot(2)+=f2(2,j);
		ftot(3)+=f2(3,j);
	}
	
	ftot1=0;
	if (ftot(1)==0 || ftot(2)==0 || ftot(3)==0)
		cerr << "\n\n ftot=0!!!" << endl;
	for (j=1; j<=N2; j++)
	{		
		f2(1,j)/=ftot(1);
		f2(2,j)/=ftot(2);
		f2(3,j)/=ftot(3);
		ftot1+=f2(1,j);
	}
	
	if (flagbug == 1)
	{
		debug << "\n p2:\n\n" << p2 << endl;
		debug << "\n l2:\n\n" << l2 << endl;
		debug << "\n Ke:\n\n" << Ke << endl;
		debug << "\n p2m: " << p2m << endl;
		debug << "\n l2m: " << l2m << endl;
		debug << "\n E2m: " << Em2 << endl;
		debug << "\n E2:\n\n" << E2 << endl;
		debug << "\n f2:\n\n" << f2 << endl;
	}
}

//--------------------------------------------------------------------------------------------
 
void PROFIT(void)			// Calcolo domanda di lavoro, profitti e net worth imprese

{
	 //cout << "\nProfit..." <<endl;
		
	 Pitot1=0;
	 Wtot1=0;
	 Pitot2=0;
	 Wtot2=0;
	 cpi(1)=0;
	 c2tot=0;
	 flag=0;
	 Qtot2=0;
	 Divtot(1)=0;
	 div1=0;
	 div2=0;
	 dNmtot=0;
	 Tax=0;
	 CF=0;
	 tax1=0;
	 	 
	 for (i=1; i<=N1; i++)									// Industria macchinari
	 {
	     Pi1(i)=S1(1,i)-c1(i)*Q1(i)-RD(2,i)+r*W1(2,i)-r*Deb2(2,i);
		if (Pi1(i) > 0)
		 {
			 div1(i)=d1*Pi1(i);
			 Divtot(1)+=div1(i);
			 if (W1(1,i) <= 0)
				 cerr << "\n\n ERRORE: W1 < 0 Prima tasse" << W1(1,i);
			 tax1=aliq*Pi1(i);
			 W1(1,i)-=tax1;
			 if (W1(1,i) <= 0 && flagdieW==0)					// imprese posso morire x motivi finanziari
			 {
				 cerr << "\n\n ERRORE flagdieW!!! " << endl;
			 }
			 if (W1(1,i) <= 0 && flagdieW==1)					// imprese non possono morire x motivi finanziari
			 {
				 W1(1,i)+=tax1;
				 tax1=0;				 
			 }
			Tax+=tax1;											// tasse totali	*/
		 }	
		 W1(1,i)=W1(1,i)+r*W1(2,i)-r*Deb2(2,i)-div1(i);		// NOTA: S1 e RD tolti prima!
		 
		 if (W1(1,i) > 1 && W1(1,i) > W1(2,i) && Deb1(1,i) > 0)
		 {
			 if (Deb1(1,i) < W1(1,i) - W1(2,i))
			 {
				 W1(1,i)-=Deb1(1,i);
				 Deb1(1,i)=0;
			 }
			 else 
			 {
				 Deb1(1,i)=Deb1(1,i) - (W1(1,i) - W1(2,i));
				 W1(1,i)=W1(2,i);
			 }
		 }
		 
		 if (flagdieW == 0)  // Se = 0, impresa pu‚àö‚â§ morire per motivi finanziari
		 {
			 if (W1(1,i) <= 0)
				 cout << "\n Morte x motivi finanziari!!!" << endl;
			 if (W1(1,i) <= 0 && debres1(i) + W1(1,i) > 0)		   
			 {
				 Deb1(1,i)-=W1(1,i);
				 debres1(i)=0;
				 W1(1,i)=1;
			 }
			 if (W1(1,i) <= 0 && debres1(i) + W1(1,i) <= 0)
			 {											   
				 Deb1(1,i)-=W1(1,i);
				 debres1(i)=0;
				 W1(1,i)=0;
			 }
		 }
		 
		 if (W1(1,i) <= 0)
		 {
			 cerr << "\n\n ERRORE: W1 < 0 Dopo tasse " << W1(1,i);
			 W1(1,i)+=0.01;
		 }

		 Pitot1+=Pi1(i);								// Profitti industria 1
		 Wtot1+=W1(1,i);								// Ricchezza industria 1
		 //size1(1,i)=Ld1(i);							// Size imprese 1
		 size1(1,i)=Q1(i);
	 }

	 if (flagbug == 1)
	 {
		 debug << "\n Q1: \n\n" << Q1 << endl;
		 debug << "\n c1: \n\n" << c1 << endl;
		 debug << "\n p1: \n\n" << p1 << endl;
		 debug << "\n S1: \n\n" << S1 << endl;
		 debug << "\n RD: \n\n" << RD << endl;
		 debug << "\n Pi1:\n\n" << Pi1 << endl;
		 debug << "\n W1:\n\n" << W1 << endl;
	 }

	 if (flagC==0)
		 G(1)=0;
	 else 
	 {
		 if (flagC==1)
	 		 G(1)=(1+0.005)*G(2);
	 	 else 
		 {
			 if (flagC==2 && LS>LD)
	 	 		 G(1)=(LS-LD)*(w(2)*wu);
	 		 else G(1)=0;
			 Gtot+=G(1);
		 }
	 }
	 	
	 if (flagCN==1)
		 Cons=w(2)*LD+Divtot(2)+G(1)+Cpast*(1+r);
	 else
		 Cons=w(2)*LD+Divtot(2)+G(1);
	 
	 if (flagTAX==1)			// Tasse su consumi
	 {
		 Cons-=aliq*w(2)*LD;	// consumi netti
		 Tax+=aliq*w(2)*LD;		// tasso consumi
	 }

	 for (j=1; j<=N2; j++)								// Industria manifatturiera
	 {
		 if (Q2(1,j) > K(j))
			 cerr << "\n ERRORE: Produzione > Capitale!" << endl;
		 cpi(1)+=p2(1,j)*f2(1,j);						// Consumer price index

		 c2tot+=c2e(j)*f2(1,j);							// Costi totali
		 if (p2(1,j)<=0 && Ke(j)==0)
			 cerr << "\n ERRORE: p2 = 0!" << endl;
	 }
	if (cpi(1) < 0.01)										// cpi >=1 per non fare
	{
		cerr << "\n\n ERRORE CPI < 0.01!!!  " << cpi(1) << endl; // esploredere variabili reali
	}

	 ALLOC();											// Alloco consumi ad imprese

	 Ntot=0;
	 dNtot=0;
	 dNmtot=0;
	
	 for (j=1; j<=N2; j++)
	 {
		 if (flagENTRY >= 2 && Ke(j) > 0)
			 ueff(j)=0;
		 else
			 ueff(j)=Q2(1,j)/K(j);							// % utilizzo capacit‚àö‚Ä† produttiva
		 Qtot2+=Q2(1,j);
		 //size2(1,j)=Ld2(j);
		 if (D2(1,j) >= Q2(1,j) + N(2,j))
		 {
			 S2(1,j)=p2(1,j)*(Q2(1,j)+N(2,j));
			 size2(1,j)=Q2(1,j)+N(2,j);
		 }												
		 else 
		 {
			 S2(1,j)=p2(1,j)*D2(1,j);
			 size2(1,j)=D2(1,j);
		 }
		 
		 if (Ke(j) > 0 && size2(1,j)>0)
			 cerr <<"\n ERRORE: Ke>0 e Size>0!!!" << endl;
		 
		 mol(j)=S2(1,j)-c2e(j)*Q2(1,j);
		 Pi2(j)=mol(j)-r*Deb2(2,j)+r*W2(2,j);
		 if (Pi2(j) > 0)
		 {
			 div2(j)=d2*Pi2(j);
			 Divtot(1)+=div2(j);
			 Tax+=aliq*Pi2(j);							// tasse
			 CF(j)-=aliq*Pi2(j);
		 }
		 
		 N(1,j)=Q2(1,j)+N(2,j)-D2(1,j);					// Scorte impresa j
		 Nm(1,j)=N(1,j)*p2(1,j);
		 if (N(1,j)<0)
		 {
			 N(1,j)=0;
			 Nm(1,j)=0;
		 }
		 Ntot+=N(1,j);									// Scorte totali
		 dN(j)=N(1,j)-N(2,j);
		 dNm(j)=Nm(1,j)-Nm(2,j);
		 dNtot+=dN(j);
		 dNmtot+=dNm(j);
		 
		 CF(j)=Pi2(j)+c2e(j)*Q2(1,j)-div2(j);

		 if (Deb2(1,j) < 0)
			 cerr << "\n\n ERRORE: debito negativo!" << endl;

		 if (CF(j) < 0)										// Cash Flow Negativo:
		 {
			 if (W2(1,j) >= -CF(j))							// Riduco ricchezza
				 W2(1,j)+=CF(j);
			 else										
			 {
				 if (Debmax2(j)-Deb2(1,j) >= -CF(j)-W2(1,j))	// Mi indebito
				 {
				 	 Deb2(1,j)-=CF(j)+W2(1,j);
				 	 W2(1,j)=1;
				 }
			 	 else										// Muoio
				 	 W2(1,j)+=CF(j);
			 }
		 }
		 else 
		 {
			 if (Deb2(1,j) == 0)
				 W2(1,j)+=CF(j);
			 else 
			 {
				 if (CF(j) > Deb2(1,j))
				 {
					 W2(1,j)+=CF(j)-Deb2(1,j);
					 Deb2(1,j)=0;
				 }
				 else
					 Deb2(1,j)-=CF(j);
			 }
		 }

		 Pitot2+=Pi2(j);								// Profitti industria 2
		 Wtot2+=W2(1,j);								// Ricchezza industria 2
		 
	 }
	if (echo==1 && (Wtot2 < 20000 || Wtot1 < 5000))
		cerr << "\n\n Ricchezza netta in esaurimento!" << endl;

	if (flagbug == 1)
	{
		debug << "\n Q2:\n\n" << Q2 << endl;
		debug << "\n D2:\n\n" << D2 << endl;
		debug << "\n N2:\n\n" << N << endl;
		debug << "\n S2:\n\n" << S2 << endl;
		debug << "\n mol:\n\n" << mol << endl;
		debug << "\n Pi2:\n\n" << Pi2 << endl;
		debug << "\n CF:\n\n" << CF << endl;
		debug << "\n W2:\n\n" << W2 << endl;
		debug << "\n Deb:\n\n" << Deb2 << endl;
	}
}

//--------------------------------------------------------------------------------------------

void ALLOC(void)								// Alloca domanda beni di consumo
{
	n=1;
	ftot=0;
	if (cpi(1)>0)
		Cres=Cons/cpi(1);
	else 
	{
		cout << "\n\n ERRORE: cpi=0!!!     " << cpi(1) << endl;
	}
	D2tot=0;
	Cpast=0;
	
	for (j=1; j<=N2; j++)
	{
		D2(1,j)=0;
		f_temp2(j)=f2(1,j);
		ftot(1)+=f_temp2(j);
		Q2temp(j)=Q2(1,j)+N(2,j);
	}

	while (Cres >= 1 && ftot(1) > 0)			// Cicla finch‚àö¬Æ tutta la domanda non ‚àö¬Æ assegnata
	{											// (Cres=0) o finch‚àö¬Æ tutte le imprese razionano 	
		Cresbis=Cres;							// (se net worth insufficiente, ftot=0)
		for (j=1; j<=N2; j++)				
		{
			if (f_temp2(j) > 0)
			{
				D_temp2(j)=Cres*f_temp2(j);
				if (n==1)
				{
					D2(1,j)+=D_temp2(j);
				}
				if (D_temp2(j) <= Q2temp(j))	// Imprese soddisfano tutta la domanda
				{
					if (n > 1)
						D2(1,j)+=D_temp2(j);
					Cresbis-=D_temp2(j);
					Q2temp(j)-=D_temp2(j);
					if (n==1) 
						l2(j)=1;				// Domanda non soddisfatta assegnata 1‚Äö√†√ª giro while
				}								// l1=1 per via dei logaritmi
				else							// Imprese razionano
				{
					if (n > 1)
						D2(1,j)+=Q2temp(j);
					Cresbis-=Q2temp(j);
					f_temp2(j)=0;			// assegnazione domanda residua
					if (n==1) 
						l2(j)=1+D_temp2(j)-Q2temp(j);
					Q2temp(j)=0;
				}
			}
		}
		ftot(1)=f_temp2.Sum();	
	    f_temp2/=ftot(1);
		Cres=Cresbis;
		n++;
	}
	Cpast=Cres;
}

//--------------------------------------------------------------------------------------------

void MACRO(void)						// Calcolo variabili macro
 
 {
	 //cout << "\n\n Macro..." <<endl;
	 	 
	 Debtot2=0;
	 Debtot1=0;
	 Am(1)=0;
	 EItot=0;
	 SItot=0;
	 Ir=0; 
	 Imon=0;
	 Cmon=0;
	 ppi(1)=0;
	 Am1=0;
	 Am2=0;
	 Am1bis=0;
	 Q2tot=0;
	 D2tot=0;
	 Netot=0;
	 Nt2tot=0;
	 Nt1tot=0;
	 A_mi=0;
	 A1_mi=0;
	 A_sd=0;
	 H1=0;
	 H2=0;
	 HI=0;
	 HIn=0;
	 DF1=0;
	 DF2=0;
	 IxHI=0;
	 IxHIn=0;
	 
	 Qtot1=Q1.Sum();
	 if (Qtot1==0)
		 contQtot1++;
	 N2eff=0;
	 N2lI=0;
	 N2hI=0;
	 N2lEI=0;
	 N2hEI=0;

	 for (j=1; j<=N2; j++)
	 {
		 Q2tot+=Q2(1,j);
		 D2tot+=D2(1,j);
		 Netot+=Ne(j);
		 Nt2tot+=N(2,j);
		 Nt1tot+=N(1,j);		 
		 Debtot2+=Deb2(1,j);
		 EItot+=EI(1,j);
		 SItot+=SI(1,j);
		 Am2+=A2e(j)*f2(1,j);
		 if (LD2>0)
		 {
			 Am(1)+=Ld2(j)/LD2*A2e(j);						// Dividere per LD per conteggiare anche lavoratori R&D
		 }													// Dividere per LD2 per non conteggiare lavoratori R&D
		 //else cerr << "\n\n ERRORE LD2==0!!!" << endl;
		 DF2+=fabs(f2(1,j)-f2(2,j));							// Calcolo indice variabilit‚àö‚Ä† quote mercato
		 if (flagENTRY< 2)
		 {
			 A_mi+=log(A2(j));
			 H2+=f2(1,j)*f2(1,j);							// Herfindahl index
		 }
		 else 
		 {
			 if (Ke(j)==0)
			 {
				 A_mi+=log(A2(j));
				 H2+=f2(1,j)*f2(1,j);						// Herfindahl index
				 N2eff++;
				 if (I(1,j)/K(j)<0.02)
					 N2lI++;
				 if (I(1,j)/K(j)>0.35)
					 N2hI++;
				 if (EI(1,j)/K(j)<0.02)
					 N2lEI++;
				 if (EI(1,j)/K(j)>0.35)
					 N2hEI++;
			 }
		 }
		 Cmon+=p2(1,j)*Q2(1,j);
	 }
	 
	 if (flagENTRY<2)
	 {
		 A_mi/=N2;		
		 H2=(H2-1/N2)/(1-1/N2);								// Riscalo Herfindahl
	 }
	 else 
	 {
		 A_mi/=N2eff;
		 H2=(H2-1/N2eff)/(1-1/N2eff);						// Riscalo Herfindahl
		 N2lI/=N2eff;
		 N2hI/=N2eff;
		 N2lEI/=N2eff;
		 N2hEI/=N2eff;
	 }
	 
	 EItot/=dim_mach;
	 SItot/=dim_mach;

	 for (j=1; j<=N2; j++)
	 {
		 if (flagENTRY < 2)
		 {
			 A_sd+=(log(A2(j))-A_mi)*(log(A2(j))-A_mi);
		 }
		 else
		 {
			 if (Ke(j) == 0)
			 {
				 A_sd+=(log(A2(j))-A_mi)*(log(A2(j))-A_mi);
			 }
			 
		 }
	 }
	 
	 if (flagENTRY < 2)
		 A_sd=sqrt(A_sd/N2);								// Varianza produttività
	 else
		 A_sd=sqrt(A_sd/N2eff);

	 for (i=1; i <=N1; i++)
	 {
		 if (Qtot1 > 0)
			 f1(1,i)=Q1(i)/Qtot1;
		 else f1(1,i)=f1(2,i);
		 H1+=f1(1,i)*f1(1,i);								// Herfindahl
		 DF1+=fabs(f1(1,i)-f1(2,i));						// Indice variabilità quote mercato
	     A1_mi+=log(A1p(i));
		 Debtot1+=Deb1(1,i);
		 Imon+=p1(1,i)*Q1(i);
		 if (Qtot1 > 0)
		 {
			 ppi(1)+=p1(1,i)*Q1(i)/Qtot1;
			 Am1+=A1p(i)*Q1(i)/Qtot1;
			 Am1bis+=A1(i)*Q1(i)/Qtot1;
		 }
		 else 
		 {
			 //cout << "\n\n Qtot1==0!" << endl;
			 ppi(1)+=p1(1,i);
			 Am1+=A1p(i);
			 Am1bis+=A1(i);
		 }
		 if (LD2>0)
		 {
			 Am(1)+=Ld1(i)/LD2*A1p(i);					// Dividere per LD per conteggiare anche lavoratori R&D
		 }												// Dividere per LD2 per non conteggiare lavoratori R&D
		 //else cerr << "\n\n ERRORE: LD2=0!!!" << endl;
	 }
	 
	 
	 if (flagFRONT==1 && flagFRONTEX==3)
	 {
		 A1avg=0;
		 A1pavg=0;
		 A1favg=0;
		 A1pfavg=0;
		 A1fmax=0;
		 A1pfmax=0;
		 for (i=1; i<=N1; i++)
		 {
			 A1avg+=A1(i)/N1;
			 A1pavg+=A1p(i)/N1;
		 }
		 for (i=1; i<=N1f; i++)
		 {
			 if (A1f(i)>A1fmax)
				 A1fmax=A1f(i);
			 if (A1pf(i)>A1pfmax)
				 A1pfmax=A1pf(i);
			 A1favg+=A1f(i)/N1f;
			 A1pfavg+=A1pf(i)/N1f;
		 }
		 A1fmin=A1fmax;
		 A1pfmin=A1pfmax;
		 A1fstd=0;
		 A1pfstd=0;
		 for (i=1; i<=N1f; i++)
		 {
			 if (A1f(i)<A1fmin)
				 A1fmin=A1f(i);
			 if (A1pf(i)<A1pfmin)
				 A1pfmin=A1pf(i);
			 A1fstd+=(A1f(i)-A1favg)*(A1f(i)-A1favg)/N1f;
			 A1pfstd+=(A1pf(i)-A1pfavg)*(A1pf(i)-A1pfavg)/N1f;
		 }
		 A1fstd=sqrt(A1fstd);
		 A1pfstd=sqrt(A1pfstd);
	 }
	 
	 if (ppi(1) < 0.01)
		 cout << "\n\n Errore: ppi < 0.01!!! " << ppi(1) << endl;
	 
	 A1_mi/=N1;
	 H1=(H1-1/N1)/(1-1/N1);									// Riscalo Herfindahl

	 if (Am(1) <= 0)
		 Am(1)=Am(2);
	 
	 WRITEPROD();										// Scrive file produttivit‚àö‚Ä†
	 
	 Def=G(1)-Tax;
	 Deb+=Def;

	 Ir=Qtot1;
	 	
	// for (j=1; j<=N2; j++)									// Calcolo Herfindahl investimenti
//	 {
//		 if (Ke(j) == 0)
//		 {
//			 IxHI+=I(1,j);
//			 IxHIn+=EI(1,j);
//		 }
//		 HI+=((I(1,j)/dim_mach)/Ir)*((I(1,j)/dim_mach)/Ir);
//		 HIn+=((EI(1,j)/dim_mach)/EItot)*((EI(1,j)/dim_mach)/EItot);
//	 }
//	 for (j=1; j<=N2; j++)									// Calcolo Herfindahl investimenti
//	 {
//		 if (Ke(j) == 0)
//		 {
//			 HI+=(I(1,j)/IxHI)*(I(1,j)/IxHI);
//			 HIn+=(EI(1,j)/IxHIn)*(EI(1,j)/IxHIn);
//		 }
//	 }
//	 if (Ir==0 || IxHI==0)
//		 HI=0;
//	 if (EItot==0 || IxHIn==0)
//		 HIn=0;
	 	 
	 if (cpi(1)>0)
	 {
		 Creal=Cons/cpi(1);
	 }
	 else
	 {
		 cerr << "\n\n ERRORE cpi==0!!!     " << cpi(1) << endl;
	 }
	 
	 if (Cpast/Creal >= 0.5)
	 {
		 contCC++;
	 }
	 
	 //cout << "\n\n Cpast: " << Cpast << endl;
	 
	 if (flagCN==1)
		 Creal-=Cpast;
	 Y=Creal+Ir;
	 GDP(1)=Creal+Ir+dNtot;
	 if (GDP(1) < 0 && GDP(1) > -1)
		 GDP(1)=0;
	 if (LD>0)
		 TFP=GDP(1)/LD;
	 else TFP=GDP(1);
	 
	 GDP1=Qtot1+Qtot2;
	 GDPm=Cmon+Imon+dNmtot;
	 
	 if (flagGDP==1)
	 {
	 GDP1=Qtot1*Imon/GDPm+Qtot2*(Cmon+dNmtot)/GDPm;					// calcolo PIL reale alla GD, dando peso a I
	 Creal=Creal*Cmon/GDPm;
	 Ir=Ir*Imon/GDPm;
	 dNtot=dNtot*dNmtot/GDPm;
	 GDP(1)=Creal+Ir+dNtot;
	 }

 	 Pitot=Pitot1+Pitot2;									// Profitti totali
 	 Wtot=Wtot1+Wtot2;										// Ricchezza netta totale
	 U(1)=(LS-LD)/LS;										// Disoccupazione 
	 if (U(1)<0)
	 {
		 umin=U(1);
		 if (U(1) < -0.00001)
			 cerr << "\n\n U < 0!!!     " << U(1) << endl;			 
	 }
	 if (t>100)
		 Umean+=U(1);
	 if (t==1)
		 cout << "\n\n Disoccupazione: " << U(1) << endl;
	 if (c2tot>0)
		 Mutot=(cpi(1)-c2tot)/c2tot;							// Mark-up globale
	 else cerr << "\n\n ERRORE c2tot=0!!!" << endl;
	 GDPonPi2=Pitot2/GDP(1);
	 GDPonPi1=Pitot1/GDP(1);
	 GDPonPi=Pitot/GDP(1);
	 dNonGDP=dNtot/GDP(1);
	 IonK=(Pitot2-(Ir2*ppi(1)))/(Ir2*ppi(1));
	 if (GDPm>1)											// Rapporto tra GDP
	 {														// nominale e deficit
		 DefonGDP=Def/GDPm;									// debito
		 DebonGDP=Deb/GDPm;
		 GonGDP=Gtot/GDPm;
	 }
	 else 
	 {
		 DefonGDP=Def;
		 DebonGDP=Deb;
		 GonGDP=Gtot;
	 }

	 if (flagTC==0 || flagTC==2 || t==1)
	 {
		 cout << "\n Itot: " << Ir << endl;
		 cout << " D2tot: " << D2tot << endl;
		 cout << " Netot: " << Netot << endl;
		 cout << " N_t-1: " << Nt2tot << endl;
		 cout << " Q2tot: " << Q2tot << endl;
		 cout << " N_t: " << Nt1tot << endl;
		 cout << " LS0: " << LS0 << endl;
		 cout << " LD2: " << LD2tot << endl;
		 cout << " LD1rd: " << LD1rdtot << endl;
		 cout << " LD1: " << LD1tot << endl;
		 //cout << " Pi2 / GDP: " << GDPonPi2 << endl;
		 //cout << " Pi1 / GDP: " << GDPonPi1 << endl;
		 cout << " Pi / GDP: " << GDPonPi << endl;
		 cout << " dN / GDP: " << dNonGDP << endl;
		 cout << " I / K : " << IonK << endl;
		 cout << " ppi: " << ppi(1) << endl;
	 }

	 if (U(2) < ustar)
		U(2)=ustar;
	 d_U=(U(1)-U(2))/U(2);
	 if (cpi(2)!=0)
		 d_cpi=(cpi(1)-cpi(2))/cpi(2);
	 else 
	 {
		 cerr << "\n\n ERRORE cpi2=0!!!     " << cpi(2) << endl;
	 }
	 if (Am(2)!=0)
		 d_Am=(Am(1)-Am(2))/Am(2);
	 else cerr << "\n\n ERRORE Am=0!!!" << endl;
	
	 if (d_cpi > cpi(2)*mdw)				// Limito variabilità componenti salario ed evito
		 d_cpi=cpi(2)*mdw;					// divisioni per zero!
	 if (d_Am > Am(2)*mdw)
		 d_Am=Am(2)*mdw;
	 if (d_U > U(2)*mdw)
		 d_U=U(2)*mdw;
	 if (d_cpi < -cpi(2)*mdw)				
		 d_cpi=-cpi(2)*mdw;					
	 if (d_Am < -Am(2)*mdw)
		 d_Am=-Am(2)*mdw;
	 if (d_U < -U(2)*mdw)
		 d_U=-U(2)*mdw;
	 
 	 dw= psi1*d_cpi + psi2*d_Am - psi3*d_U;					// Salario nominale
  	 //dw= psi1*d_cpi + psi2*d_Am - psi3*(U(1)-ustar);					// Salario nominale
	 //	 if (dw < -0.5)
//		 dw=-0.5;
//	 if (dw>0.5)
//		 dw=0.5;
	 
	 if (flagSticky==0)
		 w(1)=w(2)*(1+dw);
	 else 
	 {
		 wstar(1)=wstar(2)*(1+dw);
		 rnd=ran1(p_seed);
		 if (rnd > sticky)
			 w(1)=wstar(1);
		 else w(1)=w(2);			 
	 }

	 if (w(1) < w_min-0.001)
	 {	 
		 cerr << "\n ERRORE: salario inferiore a livello sussistenza!!!" << endl;
		 cout<< "\n w: "<<w(1) << endl;
		 w(1)=w_min;										// Salario sussistenza
	 }
	 diff_w(1)=log(w(1))-log(w(2));
	 diff_cpi(1)=log(cpi(1))-log(cpi(2));
	 
	 if (cpi(1) > 0)
		 rw=w(2)/cpi(1);										// Salario reale
	 else 
	 {
		 cerr << "\n\n ERRORE: cpi = 0!!!     " << cpi(1) << endl;
	 }
	 LS*=(1+eta);											// Offerta lavoro
	 
	 if (flagbug == 1)
	 {
		 debug << "\n Am:" << Am << endl;
		 debug << "\n A_mi:" << A_mi << endl;
		 debug << "\n A_sd:" << A_sd << endl;
		 debug << "\n w:" << w << endl;
	 }
 }

//--------------------------------------------------------------------------------------------

void WRITEPROD(void)									// Scrive i file di produttivit‚àö‚Ä†
{
	if(imc<=3)
	{
		ofstream inv_prodall1(nomefile6,ios::app);
		inv_prodall1.setf(ios::fixed);
		inv_prodall1.precision(4);
		inv_prodall1.setf(ios::right);
		if (t>1)
			inv_prodall1 << "\n";
		for (i=1; i<=N1; i++)
		{
			inv_prodall1.width(20);
			inv_prodall1 << A1p(i);
		}
		inv_prodall1.close();
		
		ofstream inv_prodall2(nomefile7,ios::app);
		inv_prodall2.setf(ios::fixed);
		inv_prodall2.precision(4);
		inv_prodall2.setf(ios::right);
		if (t>1)
			inv_prodall2 << "\n";
		for (j=1; j<=N2; j++)
		{
			inv_prodall2.width(20);
			if (Ke(j)==0)
				inv_prodall2 << A2(j);
			else inv_prodall2 << "NaN";
		}
		inv_prodall2.close();
	}
	
	if(imc<=3 && t>=101) //501)								// Scrivo produttivit‚àö‚Ä† imprese (ultimi 500 turni)
	{													// per analisi persistenza
		ofstream inv_prod2(nomefile4,ios::app);
		inv_prod2.setf(ios::fixed);
		inv_prod2.precision(4);
		inv_prod2.setf(ios::right);
		if (t>101)//(t>501)
			inv_prod2 << "\n";
		for (j=1; j<=N2; j++)
		{
			inv_prod2.width(20);
			if (f2(1,j) >= exit2 && W2(1,j) > 0 && Ke(j)==0) // flagA2(j)==0)
			{
				A2scr=log(A2(j))-A_mi;
				inv_prod2 << A2scr;
				if (Ke(j)>0)
					cerr << "\n ERRORE: Ke > 0 nel calcolo A!!!" << endl;
			}
			else 
			{
				inv_prod2 << "NaN";
				flagA2(j)=1;
			}
		 }
		inv_prod2.close();
		
		ofstream inv_prod1(nomefile5,ios::app);
		inv_prod1.setf(ios::fixed);
		inv_prod1.precision(4);
		inv_prod1.setf(ios::right);
		if (t>101)//501)
			inv_prod1 << "\n";
		for (i=1; i<=N1; i++)
		{
			inv_prod1.width(20);
			if (nclient(i) >= 1 && W1(1,i) > 0)// && flagA1(i)==0)
			{
				A1scr=log(A1p(i))-A1_mi;
				inv_prod1 << A1scr;
			}
			else 
			{
				inv_prod1 << "NaN";
				flagA1(i)=1;
			}
		 }
		inv_prod1.close();
	 }	
}

//--------------------------------------------------------------------------------------------
 
void ENTRYEXIT(void)									// Entrata - uscita imprese
{
	//if (t >= 501)
	//	nsize++;
	//else nsize=0;
	
	//GENFILESIZE1("ZS1_",nsize);
	ofstream inv_size1(nomefile0);
	inv_size1.setf(ios::fixed);
	inv_size1.precision(4);
	inv_size1.setf(ios::right);
	//GENFILESIZE2("ZS2_",nsize);
	ofstream inv_size2(nomefile3);
	inv_size2.setf(ios::fixed);
	inv_size2.precision(4);
	inv_size2.setf(ios::right);
	
	//cout << "\n\n Entry Exit!!" << endl;

	flag=2;
	ee1=0;
	ee2=0;
	W1m=0;
	W2m=0;
	Km=0;
	nwm1=0;
	nwm2=0;
	next1=0;
	next2=0;
	ind_ee1=0;
	ind_ee2=0;
	die=0;
	Kexp=0;
	f2exit=0;
	n2exit=0;	

	for (i=1; i<=N1; i++)
	{
		if (nclient(i) >= 1 && W1(1,i) > 0)
		{		
			W1m+=W1(1,i);
			nwm1++;						// Numero imprese sopravvissute industria 1
		}
		else
		{
			if (nclient(i)>=1 && W1(1,i) <= 0)
				cout << "\n Uscita per W1! " << W1(1,i) << endl;

			ind_ee1(i)=1;				// Se impresa muore ee1 = 1
			nclient(i)=0;
			W1(1,i)=0;
		}
		if (t >= 500 && size1(1,i) > 0 && size1(2,i) > 0)
		{											// Salvo size imprese con size > 0
			inv_size1.width(25);					// Nota: ultimo file‚ àö¬Æ vuoto, per cui non
			inv_size1 << size1(1,i);				// ‚àö¬Æ considerato da Matlab
			inv_size1.width(25);
			inv_size1 << size1(2,i);
			inv_size1 << "\n";
		}
	} 
	if (nwm1 > 0)
		W1m/=nwm1;
	else cerr << "\n\n ERRORE nwm1=0!" << endl;
	
	for (j=1; j<=N2; j++)
	{	
		if (f2(1,j) >= exit2 && W2(1,j) > 0)
		{
			Km+=K(j);
			W2m+=W2(1,j);
			nwm2++;						// Numero imprese sopravvissute industria 2
		}
		else
		{
			if (Ke(j)==0)				// Impresa non ‚àö¬Æ entrata nel turno precedente
				ind_ee2(j)=1;			// Se impresa muore, ee2 = 1
		}
		

		if (t >= 500 && size2(1,j) > 0 && size2(2,j) > 0)
		{											// Salvo size imprese con size > 0
			inv_size2.width(25);
			inv_size2 << size2(1,j);
			inv_size2.width(25);
			inv_size2 << size2(2,j);
			inv_size2 << "\n";
		}
	}
	
	if (nwm2 > 0)
	{
		W2m/=nwm2;
		Km/=nwm2;
	}
	else cerr << "\n\n ERRORE nwm2=0!!!" << endl;
	
	if (flagENTRY>=2)								// Se entry settore 2 non ‚àö¬Æ random copy
	{												// chi entra al tempo t, dovr‚àö‚Ä† avere
		for (j=1; j<=N2; j++)						// quota mercato e domanda attesa al tempo t+1
		{											// Per cui alla fine del tempo t, diamo
			if (Ke(j)>0)							// quota di mercato proporzionale al K dei 
			{										// sopravvissuti e mettiamo un flag per far s‚àö¬®
				f2(1,j)=Ke(j)/(Km*nwm2);			// che al tempo t+1 la domanda attesa coincida
				if (f2(1,j) < exit2)				// con capacit‚àö‚Ä† produttiva desiderata
					f2(1,j)=exit2;
				f2(2,j)=f2(1,j);					
				f2(3,j)=f2(1,j);
				Kexp(j)=1;
				f2exit+=f2(1,j);
				n2exit++;
			}
		}
		if (f2exit > f2max)
			f2max=f2exit;
	}
			
	if (flagbug == 1)
	{
		debug << "\n Entry & Exit..." << endl;
		debug << "\n Ke:\n\n" << Ke << endl;
		debug << "\n Kexp:\n\n" << Kexp << endl;
		debug << "\n f2:\n\n" << f2 << endl;
		debug << "\n ee1:\n\n" << ind_ee1 << endl;
		debug << "\n ee2:\n\n" << ind_ee2 << endl;
	}
	
	Ke=0;
	
	if (nwm1 > 0)
	{
		for (i=1; i<=N1; i++)
		{
			if (nclient(i) < 1 || W1(1,i) <= 0)		// Condizione di uscita
			{
				flag=0;
				next1++;
				for (j=1; j<=N2; j++)				// Se impresa muore, cancello i suoi clienti
				{									// e la elimino dai fornitori delle imprese
					Match(j,i)=0;					// che producono beni di consumo
				}
				while (flag == 0)					// Entrata: copia di un'impresa esistente
				{
					rni=int(ran1(p_seed)*10000)%N1+1;
					if (nclient(rni) >= 1 && W1(1,rni) > 0)// && W1(1,rni) <= W1m) 
					{	
						ee1(i)=rni;
						flag=1;
					}
				}
			}
		}
		
		for (j=1; j<=N2; j++)						// Se fornitore impresa beni di consumo muore, l'impresa
		{											// rimane senza fornitore e le macchine ordinate non
				indforn=int(fornit(j));				// vengono consegnate
				if (ind_ee1(indforn) == 1)
				{
					for (i=1; i<=N1; i++)
					{
						for (tt=t0; tt<=t; tt++)
						{
							gtemp[tt-1][i-1][j-1]=g[tt-1][i-1][j-1];
												}
					}
					if (flagENTRY>=2 && Kexp(j)==1)	// Se impresa era appena entrata, rimane senza K...
					{								// ... muore
						Kexp(j)=0;
						//Ke(j)=EI(1,j);
						f2(1,j)=0;
						f2(2,j)=0;
						f2(3,j)=0;
						cout << "\n\n Impresa appena entrata non ha K x morte fornitore!!!" << endl;
					}
					EI(1,j)=0;
					EI(2,j)=0;
					SI(1,j)=0;
					SI(2,j)=0;
					I(1,j)=0;
					I(2,j)=0;
					fornit(j)=0;						// Fornitore assegnato turno successivo					
				}
		}
		
		for (i=1; i<=N1; i++)
		{
			if (ee1(i) > 0)
			{
				iii=int(ee1(i));
				size1(1,i)=0;
				size1(2,i)=0;
				f1(1,i)=0;
				f1(2,i)=0;
				if (flagENTRY==0 || flagENTRY==2)			// Random Copy
				{
					W1(1,i)=W1(1,iii);
					W1(2,i)=W1(2,iii);
					A1(i)=A1(iii);
					c1(i)=c1(iii);
					p1(1,i)=p1(1,iii);
					p1(2,i)=p1(2,iii);
					A1p(i)=A1p(iii);
					A(t,i)=A1(iii);
				}
				else if (flagENTRY==1 || flagENTRY==3)		// Entra impresa medio piccola
				{
					rnd=ran1(p_seed);						
					rnd=w1inf+rnd*(w1sup-w1inf);			// Estraggo da U(w1inf,w1sup)
					W1(1,i)=rnd*W1m;						// Impresa entra con frazione W media
					W1(2,i)=rnd*W1m;
					if (flagFRONT == 1)						// Frontiera Esogena
					{
						if (flagFRONTEX == 0 || flagFRONTEX == 2) // Prod. agganciata frontiera
						{
							rnd=betadev(b_a2,b_b2,p_seed);		// Estraggo da Beta(b_a,b_b)
							A1(i)=rnd*A1front(2);				// Impresa entra con frazione
							rnd=betadev(b_a2,b_b2,p_seed);		// produttività frontiera
							A1p(i)=rnd*A1pfront(2);
						}
						if (flagFRONTEX == 1 || flagFRONTEX == 3) // Prod. agganciata a incumbent
						{										// come frontiera endogena
							rnd=betadev(b_a2,b_b2,p_seed);		// Estraggo da Beta(b_a,b_b)
							A1(i)=rnd*A1top(2)*(1+uu7);			// Impresa entra con frazione
							rnd=betadev(b_a2,b_b2,p_seed);		// produttività data da
							A1p(i)=rnd*A1ptop(2)*(1+uu7);       // impresa miglore * (1+uu7)
						}
						
					}
					else if (flagFRONT == 2)				// Frontiera Endogena
					{
						rnd=betadev(b_a2,b_b2,p_seed);		// Estraggo da Beta(b_a,b_b)
						A1(i)=rnd*A1top(2)*(1+uu7);			// Impresa entra con frazione
						rnd=betadev(b_a2,b_b2,p_seed);		// produttività data da
						A1p(i)=rnd*A1ptop(2)*(1+uu7);       // impresa miglore * (1+uu7)
					}
				
					//else cerr << "\n\n ERRORE Frontiera!!!" << endl;
					A(t,i)=A1(i);
					if (A1p(i) > 0)
						c1(i)=w(2)/(A1p(i)*a);
					else cerr << "\n\n ERRORE: A1p=0!!!" << endl;
					p1(1,i)=(1+mi1)*c1(i);
					p1(2,i)=p1(1,i);
					if (A1(i) > A1max && A1p(i) > A1pmax)
					{
						contA1entry++;
						//cout << "\n Produttivita entry > incumbent!" << endl;
					}
				}
				else 
				{
					cerr << "\n\n ERRORE flagENTRY!!!" << endl;
				}
				S1(1,i)=p1(1,i)*step;//S1(1,iii);			// Se progr. tecn ‚àö¬Æ endogeno, S1 determina quanto R&D
				S1(2,i)=S1(1,i);							// fa impresa appena entrata
				stepbis=step;								// Step=N2/N1;
				while (stepbis > 0)							// Creazione lista clienti
				{											// # clienti = N2/N1
					rni=int(ran1(p_seed)*10000)%N2+1;
					if (Match(rni,i) == 0)
					{
						Match(rni,i)=1;
						stepbis--;
					}
				}
			}
		} 
	}
	else 
	{
		cerr << "\n ERRORE: nessun sopravvissuto nell'industria 1!" << endl;
	}
	
	if (next1 > nextmax1)
		nextmax1=next1;

	//flag=2;
	if (nwm2 > 0)
	{
		for (j=1; j<=N2; j++)
		{
			if (f2(1,j) < exit2 || W2(1,j) <= 0)			// Condizione di uscita
			{
				if (Kexp(j)==1) {
					cerr <<"\n ERRORE: uscita con Kexp=1!!!" << endl;
				}
				
				flag=0;
				next2++;
				indforn=int(fornit(j));
				
				if (indforn>=1)
				{
					if (flagdie==1)
						die(indforn)+=p1(1,indforn)*I(1,j);
					Match(j,fornit(j))=0;				// Impresa macchinari perde cliente
					fornit(j)=0;
				}
				while (flag == 0)						// Entrata: copia di un'impresa esistente
				{
					rni=int(ran1(p_seed)*10000)%N2+1;
					if (f2(1,rni) >= exit2 && W2(1,rni) > 0)// && W2(1,rni) <= W2m) 
					{
						ee2(j)=rni;
						flag=1;
					}
				}
			}
		}
		
		for (j=1; j<=N2; j++)
		{
			if (ee2(j) > 0)
			{
				jjj=int(ee2(j));						// Indice impresa da copiare
				N(1,j)=0;
				N(2,j)=0;
				Nm(1,j)=0;
				Nm(2,j)=0;
				size2(1,j)=0;
				size2(2,j)=0;
				EI(1,j)=0;								// Chi entra non ha investito...
				EI(2,j)=0;
				SI(1,j)=0;
				SI(2,j)=0;
				I(1,j)=0;
				I(2,j)=0;
				Deb2(1,j)=0;							// ... e non ha debiti
				Deb2(2,j)=0;
				rni=int(ran1(p_seed)*10000)%N1+1;		// Fornitore casuale
				fornit(j)=rni;
				Match(j,rni)=1;
				if (flagENTRY==0 || flagENTRY==1)		// Entata Random Copy
				{
					f2(1,j)=f2(1,jjj);
					f2(2,j)=f2(2,jjj);
					f2(3,j)=f2(2,jjj);
					mu2(1,j)=mu2(1,jjj);
					mu2(2,j)=mu2(2,jjj);
					E2(1,j)=E2(1,jjj);
					E2(2,j)=E2(2,jjj);
					W2(1,j)=W2(1,jjj);
					W2(2,j)=W2(2,jjj);
					D2(1,j)=D2(1,jjj);
					D2(2,j)=D2(2,jjj);
					D2(3,j)=D2(3,jjj);
					D2(4,j)=D2(4,jjj);
					D2(5,j)=D2(5,jjj);
					De(1,j)=De(1,jjj);
					De(2,j)=De(2,jjj);
					De(3,j)=De(3,jjj);
					c2(1,j)=c2(1,jjj);
					c2(2,j)=c2(2,jjj);
					Q2(1,j)=Q2(1,jjj);
					Q2(2,j)=Q2(2,jjj);
					S2(1,j)=S2(1,jjj);
					S2(2,j)=S2(2,jjj);					
					l2(j)=l2(jjj);				
					K(j)=K(jjj);
					p2(1,j)=p2(1,jjj);	
					p2(2,j)=p2(2,jjj);
					S2(1,j)=S2(1,jjj);				// Serve per avere limite massimo indebitamento
					S2(2,j)=S2(2,jjj);					
					mol(j)=mol(jjj);
					n_mach(j)=n_mach(jjj);
					for (i=1; i<=N1; i++)
					{
						for (tt=t0; tt<=t; tt++)	// Non avendo investito lo stock di capitale
						{							// ‚àö¬Æ costante
							g[tt-1][i-1][j-1]=g[tt-1][i-1][jjj-1];
							gtemp[tt-1][i-1][j-1]=g[tt-1][i-1][jjj-1];
							g_c[tt-1][i-1][j-1]=g_c[tt-1][i-1][jjj-1];
							age[tt-1][i-1][j-1]=0;	// Et‚àö‚Ä† macchine ‚àö¬Æ zero
						}  
					} 
				}
				else if (flagENTRY==2 || flagENTRY==3)				// Entra impresa medio piccola
				{
					rnd=ran1(p_seed);						
					rnd=w2inf+rnd*(w2sup-w2inf);					// Estraggo da U(w2inf,w2sup)
					W2(1,j)=rnd*W2m;								// Impresa entra con frazione W media
					W2(2,j)=rnd*W2m;
					rnd=ran1(p_seed);						
					rnd=kinf+rnd*(ksup-kinf);						// Estraggo da U(Kinf,Ksup)
					Ke(j)=rnd*Km;									// Impresa entra con frazione K medio
					Ke(j)=floor(Ke(j)/dim_mach)*dim_mach;			// Arrotondo per difetto
					if (Ke(j)<=0)									// Impresa deve avere almeno 1 macchina
						Ke(j)=dim_mach;
					rni=int(fornit(j));
					if (A1(rni)<=0)
						cerr << "\n\n ERRORE: A1=0!!!" << endl;
					//if (W2(1,j) < Ke(j)/dim_mach*p1(1,rni)+(u*Ke(j)*(1+theta)*w(1)/A1(rni)))
					if (W2(1,j) < Ke(j)/dim_mach*p1(1,rni)+(Ke(j)*w(1)/A1(rni)))
					{																		// W2 deve permettere di finanziare
						//Ke(j)=W2(1,j)/(p1(1,rni)/dim_mach+(u*(1+theta)*w(1)/A1(rni)));	// K e produzione primi 2 periodi
						Ke(j)=W2(1,j)/(p1(1,rni)/dim_mach+(w(1)/A1(rni)));
						Ke(j)=floor(Ke(j)/dim_mach)*dim_mach;
						if (Ke(j) <= 0)
						{
							//cout << "\n Ke=0 dopo arrotondamento! Forzato Ke=1!" << endl;
							Ke(j)=dim_mach;
							//W2(1,j)=p1(1,rni)+u*dim_mach*(1+theta)*w(1)/A1(rni);
							W2(1,j)=p1(1,rni)+dim_mach*w(1)/A1(rni);
							W2(2,j)=W2(1,j);
						}
					}
					f2(1,j)=0;
					f2(2,j)=0;
					f2(3,j)=0;
					mu2(1,j)=mi2;
					mu2(2,j)=mi2;
					E2(1,j)=0;
					E2(2,j)=0;
					D2(1,j)=0;
					D2(2,j)=0;
					D2(3,j)=0;
					D2(4,j)=0;
					D2(5,j)=0;
					De(1,j)=0;
					De(2,j)=0;
					De(3,j)=0;
					A2(j)=0;
					c2(1,j)=0;
					c2(2,j)=0;
					Q2(1,j)=0;
					Q2(2,j)=0;
					S2(1,j)=0;
					S2(2,j)=0;					
					l2(j)=1;				
					K(j)=0;
					p2(1,j)=0;	
					p2(2,j)=0;
					mol(j)=0;
					n_mach(j)=0;
					for (i=1; i<=N1; i++)
					{
						for (tt=t0; tt<=t; tt++)	
						{							
							g[tt-1][i-1][j-1]=0;
							gtemp[tt-1][i-1][j-1]=0;
							g_c[tt-1][i-1][j-1]=0;
							age[tt-1][i-1][j-1]=0;	// Et‚àö‚Ä† macchine ‚àö¬Æ zero
						}  
					} 
					
				}
				else 
				{
					cerr << "\n\n ERRORE flagENTRY!!!" << endl;
				}
			}  
		}  
		ftot=0;										// Riscaliamo quote mercato
		for (j=1; j<=N2; j++)
		{
			ftot(1)+=f2(1,j);
			ftot(2)+=f2(2,j);
			ftot(3)+=f2(3,j);
		}
		Em2(1)=0;
		for (j=1; j<=N2; j++)
		{
			f2(1,j)/=ftot(1);
			f2(2,j)/=ftot(2);
			f2(3,j)/=ftot(3);
		}
	}
	else 
	{
		cerr << "\n ERRORE: nessun sopravvissuto nell'industria 2!" << endl;
		for (j=1; j<=N2; j++)
		{
			cout << "\n Impresa " << j << "\t\t W2: " << W2(1,j) << "\t\t f2: " << f2(1,j) << endl;
			cout << " Pi2: " << Pi2(j) << "\t\t Q2: " << Q2(1,j) << "\t\t I: " << I(1,j) << endl;
		}
	}
	
	if (next2 > nextmax2)
		nextmax2=next2;

	inv_size1.close();
	inv_size2.close();
	
	nclient=0;
	for (i=1; i<=N1; i++)
	{
		for (j=1; j<=N2; j++)
		{
			nclient(i)+=Match(j,i);
		}
	}
	for (i=1; i<=N1; i++)
	{
		if (nclient(i) == 0)
		{
			stepbis=step;
			while (stepbis > 0)						// Creazione lista clienti
			{
				rni=int(ran1(p_seed)*10000)%N2+1;
				if (Match(rni,i) == 0)
				{
					Match(rni,i)=1;
					stepbis--;
				}
			}
		}
	}
		
	for (i=1; i<=N1; i++)
	{
		for (j=1; j<=N2; j++)
		{
			nclient(i)+=Match(j,i);
		}
	}

	if (flagbug == 1)
	{
		debug << "\n Entry & Exit..." << endl;
		debug << "\n ee1:\n\n" << ind_ee1 << endl;
		debug << "\n A1:\n\n" << A1 << endl;
		debug << "\n A1p:\n\n" << A1p << endl;
		debug << "\n W1:\n\n" << W1 << endl;
		debug << "\n ee2:\n\n" << ind_ee2 << endl;
		debug << "\n W2:\n\n" << W2 << endl;
		debug << "\n Ke:\n\n" << Ke << endl;
		debug << "\n f2:\n\n" << f2 << endl;
		debug << "\n fornitore:\n\n" << fornit << endl;
		debug << "\n nclient:\n\n" << nclient << endl;
		debug << "\n Match:\n\n" << Match << endl;
	}

	for (j=1; j<=N2; j++)				
	{ 
		machtemp=n_mach(j);
	 	for (i=1; i<=N1; i++)   
		{ 
			for (tt=t0; tt<=t; tt++)
			{
				machtemp-=g[tt-1][i-1][j-1];
			}
		}
		if (machtemp != 0)
		{
			cerr << "\n\n K != g[]!!!   TRE   " << machtemp << endl; 
			cout << "\n K(): " << K(j) << endl;
			cout << "\n I(): " << I(2,j)<< endl;
			cout << "\n EI(): " << EI(2,j)<< endl;
			cout << "\n SI(): " << SI(2,j)<< endl;
		}
	}
}

//--------------------------------------------------------------------------------------------

void TECHANGEX(void)						// Progresso Tecnico settore macchinari
 { 
	
	 for (i=1; i<=N1; i++)
	 {
		 epss=100;
		 while (epss < -uu2 || epss > uu2)
		 {
			epss=gasdev(p_seed);
		 }

		 if (flagTC==2)  // Progresso tecnico deterministico
			 epss=uu;

	  	 Anew(i)=A1(i)*(1+epss);			// Produttivit‚àö‚Ä† nuovo macchinario

		 if (flagTC==1 || flagTC==2)
		 {
		 	 if (Anew(i) > A1(i))
			 {
			 	 tao(i)++;						// Incremento generazione macchinari
			 	 A1(i)=Anew(i);					// Aggiorno produttivit‚àö‚Ä† macchinario impresa i
			 	 //A(tao(i),i)=Anew(i);			// Aggiorno matrice macchine esistenti
				 if (flagTC==2 && t<T)			// Aggiorno matrice macchine esistenti
					 A(t+1,i)=Anew(i);
				 else							// NOTA: ci potrebbe essere un errore: A(t+1,i) invece che A(t,i)
					 A(t,i)=Anew(i);			
			 }
			 else
			 {
				 if (t>1)
					 A(t,i)=A(t-1,i);
			 }
		 }
		 else 
		 {
			 if (t>1)
				 A(t,i)=A(t-1,i);
		 }

		 epss=100;
		 while (epss < -uu1 || epss > uu1)		// Innovazione processo
		 {
			 epss=gasdev(p_seed);
		 }

		 if (flagTC==2)
			 epss=uu;
		 
		 if ((flagTC==1 || flagTC==2))
		 {
			 if (A1p(i)*(1+epss) > A1p(i))// && i==1)
			 {
			 	 A1p(i)*=(1+epss);					// Aggiorno produttivit‚àö‚Ä† macchinario impresa i
			 }
		 }
	 }
	 if (flagbug == 1)
	 {
		 debug << "\n Exogenous Technichal Change..." << endl;
		 debug << "\n Anew:\n\n" << Anew << endl;
		 debug << "\n A1:\n\n" << A1 << endl;
	 }
 }

//--------------------------------------------------------------------------------------------

void TECHANGEND(void) // Cambiamento Tecnico Endogeno
{
	//cout << "\n\nTechange..." <<endl;
	
	if (flagFRONT == 1)
	{
		A1front(1)=A1front(2)*(1+uu3);
		A1pfront(1)=A1pfront(2)*(1+uu4);
	}
	if (flagFRONT == 1 && flagFRONTEX==3)		// Espansione frontiera esogena
	{											// formata da imprese stilizzate
		for (i=1; i<=N1; i++)
		{
			A1w(i)=A1(i);
			A1pw(i)=A1p(i);
		}
		for (i=1; i<=N1f; i++)
		{
			
			rnd=betadev(b_a3,b_b3,p_seed);					// Estraggo da Beta(b_a3,b_b3);
			rnd=uinf+rnd*(usup-uinf);						// Riscalo su supporto (uinf,usup)
			if (rnd>0)
			{
				A1f(i)=A1f(i)*(1+rnd);
				if (flagPAT>=1)									// Se ci sono patent, inizializzo durata
					Patw(N1+i)=patdur+1;
			}
			rnd=betadev(b_a3,b_b3,p_seed);					// Estraggo da Bet(b_a3,b_b3);
			rnd=uinf+rnd*(usup-uinf);							// Riscalo su supporto (uinf,usup)	
			if (rnd>0)
			{
				A1pf(i)=A1pf(i)*(1+rnd);
				if (flagPAT>=1)									// Se ci sono patent, inizializzo durata
					Patw(N1+i)=patdur+1;
			}
			A1w(N1+i)=A1f(i);									// Creo unico vettore con imprese nazionaoi e straniere
			A1pw(N1+i)=A1pf(i);
			if (flagPAT>=1 && Pat(N1+i)>=1)
			{
				Patw(N1+i)--;
			}
		}
	}
	
	Inn=0;
	Imm=0;
	A1inn=0.00001;								// Inizializzo in modo che se
	A1pinn=0.00001;								// non innovo o imito scelgo
	A1imm=0.00001;								// macchina corrente
	A1pimm=0.00001;
	Inntot1=0;
	Inntot2=0;
	Immtot1=0;
	Immtot2=0;
	parbertot=0;
	parbertot1=0;
	A1max=0;
	A1pmax=0;
	//A1fmax=0;
	//A1pfmax=0;
	
	if (flagbug == 1)
	{
		debug << "\n Endogenous Technical Change..." << endl;
		debug << "\n A1top: " << A1top << endl;
		debug << "\n A1ptop: " << A1ptop << endl;
		debug << "\n A1: \n\n" << A1 << endl;
		debug << "\n A1p: \n\n" << A1p << endl;
		debug << "\n A:\n\n" << A << endl;
		debug << "\n C:\n\n" << C << endl;
		debug << "\n w: " << w << endl;
		debug << "\n c1: \n\n" << c1 << endl;
		debug << "\n p1: \n\n" << p1 << endl;
	}
	
	for (i=1; i<=N1; i++)
	{						
		RD(1,i)=nu*S1(1,i);							// Allocazione risorse R&D				
		if (S1(1,i)==0)								// Se nel periodo precedente clienti non hanno investito,
		{											// impresa investe in RD come nel periodo precedente
			RD(1,i)=RD(2,i);						// Ci‚àö‚â§ accade anche nel periodo t+1 alle imprese
			if (nclient(i) < 1)						// entrate al tempo t (perch‚àö¬Æ c'‚àö¬Æ time-to-build)
			{
				cerr << "\n\n ERRORE: nclinet < 1!!!" << endl;
			}
		}
		if (w(1)>0)
			Ld1rd(i)=RD(1,i)/w(1);								// Domanda lavoratori R&D
		else cerr << "\n\n ERRORE: w=0!!!" << endl;
		if (flagRD==0)											// R&D nominale (utilizza fatturato)
		{
			RDin(i)=RD(1,i)*xi;
			RDim(i)=RD(1,i)*(1-xi);
		}
		else													// R&D reale (utilizza lavoratori)
		{
			RDin(i)=Ld1rd(i)*xi;
			RDim(i)=Ld1rd(i)*(1-xi);
		}
		
		parber=1-exp(-o1*RDin(i));							// Primo stadio innovazione
		parbertot+=parber;									// Parametro Bernulli
		Inn(i)=bnldev(parber,1,p_seed);						// Estraggo da Binomiale con n=1 e p=parber
		
		parber=1-exp(-o2*RDim(i));							// Primo stadio imitazione
		parbertot1+=parber;									// Parametro Bernulli
		Imm(i)=bnldev(parber,1,p_seed);						// Estraggo da Binomiale con n=1 e p=parber
		
		if (Inn(i) == 1)									// Impresa innova
		{
			Inntot1++;
			if (flagFRONT==1)								// Frontiera Esogena
			{
				if (flagFRONTEX==0 || flagFRONTEX==3)		// Innovation local
				{
					rnd=betadev(b_a1,b_b1,p_seed);			// Estraggo da Beta(b_a1,b_b1);
					rnd=uu5+rnd*(uu6-uu5);					// Riscalo su supporto (uu5,uu6)
					A1inn(i)=A1(i)*(1+rnd);
					if (A1inn(i) > A1front(1) && flagFRONTEX==1)
					{
						A1inn(i)=A1front(1);
						contfront++;
					}
					rnd=betadev(b_a1,b_b1,p_seed);			// Estraggo da Beta(b_a1,b_b1);
					rnd=uu5+rnd*(uu6-uu5);					// Riscalo su supporto (uu5,uu6)
					A1pinn(i)=A1p(i)*(1+rnd);
					if (A1pinn(i) > A1pfront(1) && flagFRONTEX==1)
					{
						A1pinn(i)=A1pfront(1);
						contfront++;
					}
				}
				if (flagFRONTEX==1 || flagFRONTEX==2)		// Innovazione agganciata frontiera
				{
					rnd=betadev(b_a1,b_b1,p_seed);					// Estraggo da Beta(b_a1,b_b1);
					A1inn(i)=rnd*A1front(1);
					rnd=betadev(b_a1,b_b1,p_seed);					// Estraggo da Beta(b_a1,b_b1);
					A1pinn(i)=rnd*A1pfront(1);
				}				
			}
				
			if (flagFRONT==2)								// Frontiera Endogena
			{
				rnd=betadev(b_a1,b_b1,p_seed);					// Estraggo da Beta(b_a1,b_b1);
				rnd=uu5+rnd*(uu6-uu5);						// Riscalo su supporto (uu5,uu6)
				A1inn(i)=A1(i)*(1+rnd);
				rnd=betadev(b_a1,b_b1,p_seed);					// Estraggo da Bet(b_a1,b_b1);
				rnd=uu5+rnd*(uu6-uu5);							// Riscalo su supporto (uu5,uu6)	
				A1pinn(i)=A1p(i)*(1+rnd);
 			}

			if (flagFRONT==1 && flagFRONTEX < 3 && A1inn(i) > A1front(1))
			{
				cerr << "\n\n ERRORE: oltre la frontiera!!!!" << " A1: " << A1inn(i) << " A1front: " << A1front(1) << endl;
				A1inn(i)=A1front(1);
				contfront++;
			}

			if (flagFRONT==1 && flagFRONTEX < 3 && A1pinn(i) > A1pfront(1))
			{
				cerr << "\n\n ERRORE: oltre la frontiera!!!!" <<  " A1p: " << A1pinn(i) << " A1pfront: " << A1pfront(1) << endl;
				A1pinn(i)=A1pfront(1);
				contfront2++;
			}
			
			if (flagPAT==2)									// Patent breadth
			{
				for (ii=1; ii<=N1; ii++)
				{
					if (A1inn(i) >= A1(ii)*(1-patbre) && A1inn(i) <= A1(ii)*(1+patbre))
					{
						A1inn(i)=0.00001;					// Non innovo
						A1pinn(i)=0.00001;
					}
					if (A1pinn(i) >= A1p(ii)*(1-patbre) && A1pinn(i) <= A1p(ii)*(1+patbre))
					{
						A1inn(i)=0.00001;
						A1pinn(i)=0.00001;					// Non innovo
					}
				}
				if (flagFRONT==1 && flagFRONTEX==3)	// Se c'è frontiera esogena stilizzata
				{									// non si possono copiare neanche quest'ultime
					for (ii=1; ii<=N1f; ii++)
					{
						if (A1inn(i) >= A1f(ii)*(1-patbre) && A1inn(i) <= A1f(ii)*(1+patbre))
						{
							A1inn(i)=0.00001;					// Non innovo
							A1pinn(i)=0.00001;
						}
						if (A1pinn(i) >= A1pf(ii)*(1-patbre) && A1pinn(i) <= A1pf(ii)*(1+patbre))
						{
							A1inn(i)=0.00001;
							A1pinn(i)=0.00001;					// Non innovo
						}
					}
				}
			}
			if (A1pinn(i)==0 || A1inn(i)==0 || A1p(i)==0 || A1(i)==0)
				cerr << "\n\n ERRORE A1=0!!!" << endl;
		}
		
		if (Imm(i) == 1)									// Impresa imita
		{		
			if (Imm(i) == 1 && xi == 1)
				cerr << "\n\n ERRORE: imitazione non richiesta!!!" << endl;
			Immtot1++;
			Tdtot=0;
			if (flagFRONT==1 && flagFRONTEX == 3)			// Imito imprese nazionali e imprese paese sviluppato
			{												// su frontiera esogena
				for (ii=1; ii<=N1+N1f; ii++)						// Calcolo distanza tecnologica con imprese frontiera
				{
					Tdw[ii]=sqrt(((A1w(ii)-A1(i))*(A1w(ii)-A1(i))) + ((A1pw(ii)-A1p(i))*(A1pw(ii)-A1p(i))));
					if (Tdw[ii]>0)
						Tdw[ii]=1/Tdw[ii];
					else Tdw[ii]=1/0.01;	
					Tdtot+=Tdw[ii];
				}
				for (ii=1; ii<=N1+N1f; ii++)						// Normalizzo e ordino su segmento (0,1)
				{
					Tdw[ii]/=Tdtot;
					Tdw[ii]+=Tdw[ii-1];	
				}
				rnd=ran1(p_seed);								// Estraggo da uniforme tra 0 e 1
				for (ii=1; ii<=N1+N1f; ii++)						// Scelgo impresa da imitare
				{
					if (rnd <= Tdw[ii] && rnd > Tdw[ii-1] && Patw(ii)==0)
					{
						A1imm(i)=A1w(ii);
						A1pimm(i)=A1pw(ii);
						if (ii <= N1)
							contIMMnaz++;
						else
							contIMMfor++;
					}
				}	
			}
			else												// Imito solo imprese local
			{
				for (ii=1; ii<=N1; ii++)						// Calcolo distanza tecnologica con altre imprese
				{
					Td[ii]=sqrt(((A1(ii)-A1(i))*(A1(ii)-A1(i))) + ((A1p(ii)-A1p(i))*(A1p(ii)-A1p(i))));
					if (Td[ii]>0)
						Td[ii]=1/Td[ii];
					else Td[ii]=0;	
					Tdtot+=Td[ii];
				}
				for (ii=1; ii<=N1; ii++)						// Normalizzo e ordino su segmento (0,1)
				{
					Td[ii]/=Tdtot;
					Td[ii]+=Td[ii-1];	
				}
				rnd=ran1(p_seed);								// Estraggo da uniforme tra 0 e 1
				for (ii=1; ii<=N1; ii++)						// Scelgo impresa da imitare
				{
					if (rnd <= Td[ii] && rnd > Td[ii-1] && Pat(ii)==0)
					{
						A1imm(i)=A1(ii);
						A1pimm(i)=A1p(ii);
					}
				}
			}
			
			if (flagPAT==2)									// Patent breadth
			{
				for (ii=1; ii<=N1; ii++)
				{
					if (A1imm(i) >= A1(ii)*(1-patbre) && A1imm(i) <= A1(ii)*(1+patbre))
					{
						A1imm(i)=0.00001;					// Non innovo
						A1pimm(i)=0.00001;
					}
					if (A1pimm(i) >= A1p(ii)*(1-patbre) && A1pimm(i) <= A1p(ii)*(1+patbre))
					{
						A1imm(i)=0.00001;
						A1pimm(i)=0.00001;					// Non innovo
					}
				}
				if (flagFRONT==1 && flagFRONTEX==3)	// Se c'è frontiera esogena stilizzata
				{									// non si possono copiare neanche quest'ultime
					for (ii=1; ii<=N1f; ii++)
					{
						if (A1imm(i) >= A1f(ii)*(1-patbre) && A1imm(i) <= A1f(ii)*(1+patbre))
						{
							A1imm(i)=0.00001;					// Non innovo
							A1pimm(i)=0.00001;
						}
						if (A1pinn(i) >= A1pf(ii)*(1-patbre) && A1pinn(i) <= A1pf(ii)*(1+patbre))
						{
							A1imm(i)=0.00001;
							A1pimm(i)=0.00001;					// Non innovo
						}
					}
				}
			}
			
			if (A1pimm(i)==0 || A1imm(i)==0 || A1p(i)==0 || A1(i)==0)
				cerr << "\n\n ERRORE A1=0!!!" << endl;
		}
		
		if ( ((1+mi1)*w(1)/(A1pimm(i)*a))*(w(1)/A1imm(i))*b < ((1+mi1)*w(1)/(A1p(i)*a))*(w(1)/A1(i))*b )
		{													// Scelgo macchina con minore p*c*b
			A1(i)=A1imm(i);									// Aggiorno produttivit‚àö‚Ä† macchinario impresa i
			A1p(i)=A1pimm(i);								// Aggiorno produttivit‚àö‚Ä† processo impresa i
			Immtot2++;
			contIMM++;
			if (flagPAT>=1)									// Se ci sono patent, non posso brevettare
			{												// perchè sto copiando
				Pat(i)=0;
				Patw(i)=0;
			}
		}		
		
		if ( ((1+mi1)*w(1)/(A1pinn(i)*a))*(w(1)/A1inn(i))*b < ((1+mi1)*w(1)/(A1p(i)*a))*(w(1)/A1(i))*b )
		{													// Scelgo macchina con minore p*c*b				
			tao(i)++;										// Incremento generazione macchinari
			A1(i)=A1inn(i);									// Aggiorno produttivit‚àö‚Ä† macchinario impresa i
			A1p(i)=A1pinn(i);								// Aggiorno produttivit‚àö‚Ä† processo impresa i
			Inntot2++;
			contINN++;
			if (flagPAT>=1)									// Se ci sono patent, inizializzo durata
			{
				Pat(i)=patdur+1;
				Patw(i)=patdur+1;
			}
		}
		
		if (t < T)			
			A(t+1,i)=A1(i);									// Aggiorno produttività matrice macchine
		
		if (flagPAT>=1 && Pat(i)>=1)
		{
			Pat(i)--;										// Riduco durata patent
			Patw(i)--;
		}
		
		if (flagPAT>2)
			cerr << "\n\n ERRORE Patent!!!" << endl;
	}	
	
	LD1rdtot=Ld1rd.Sum();									// Domanda lavoro aggregata R&D
	parbertot/=N1r;											// Parametro medio Bernulli (tienen conto anche imprese con RD=0)
	parbertot1/=N1r;
	Inntot1/=N1r;											// Percentuale imprese innovano
	Inntot2/=N1r;											// Percentuale imprese innovano e scelgono nuova macchina
	Immtot1/=N1r;											// Percentuale imprese imitano
	Immtot2/=N1r;											// Percentuale imprese imitano macchina migliore
	
	//if (t==100 || t==200 || t==300 || t==400 || t==500 || t==600)
	//	{
	//		cout << "\n\n Par. medio Bernulli Inn: " << parbertot << endl;;
	//		cout << " Innovazioni 1 stadio: " << Inntot1 << endl;
	//		cout << " Innovazioni 2 stadio: " << Inntot2 << endl;
	//		cout << "\n\n Par. medio Bernulli Imm: " << parbertot1 << endl;;
	//		cout << " Imitazioni 1 stadio: " << Immtot1 << endl;
	//		cout << " Imitazioni 2 stadio: " << Immtot2 << endl;
	//	}
	

	A1top(1)=1;											// Calcolo top incumbent
	A1ptop(1)=1;
	for (i=1; i<=N1; i++)
	{
		if (A1(i) > A1top(1))
			A1top(1)=A1(i);
		if (A1p(i) > A1ptop(1))
			A1ptop(1)=A1p(i);
	}
	if ((A1top(1)==1 || A1ptop(1)==1) && t>1)
	{
		cerr << "\n\n ERRORE: A1top o A1ptop == 1 !!!" << endl;
	}
	
	for (i=1; i<=N1; i++)									// Calcolo max produttività imprese esistenti
	{
		if (A1(i)>A1max)
			A1max=A1(i);
		if (A1p(i)>A1pmax)
			A1pmax=A1p(i);
	}
	
//	A1favg=0;
//	A1pfavg=0;
//	if (flagFRONT==1 && flagFRONTEX==3)
//	{
//		for (i=1; i<=N1f; i++)
//		{
//			if (A1f(i)>A1fmax)
//				A1fmax=A1f(i);
//			if (A1pf(i)>A1pfmax)
//				A1pfmax=A1pf(i);
//			A1favg+=A1f(i)/N1f;
//			A1pfavg+=A1pf(i)/N1f;
//		}
//		A1fmin=A1fmax;
//		A1pfmin=A1pfmax;
//		A1fstd=0;
//		A1pfstd=0;
//		for (i=1; i<=N1f; i++)
//		{
//			if (A1f(i)<A1fmin)
//				A1fmin=A1f(i);
//			if (A1pf(i)<A1pfmin)
//				A1pfmin=A1pf(i);
//			A1fstd+=(A1f(i)-A1favg)*(A1f(i)-A1favg)/N1f;
//			A1pfstd+=(A1pf(i)-A1pfavg)*(A1pf(i)-A1pfavg)/N1f;
//		}
//		A1fstd=sqrt(A1fstd);
//		A1pfstd=sqrt(A1pfstd);
//	}
		
	if (flagbug == 1)
	{
		debug << "\n\n S1: \n\n" << S1 << endl;
		debug << "\n\n R&D: \n\n" << RD << endl;
		debug << "\n\n Ld1rd: \n\n" << Ld1rd << endl;
		debug << "\n\n Inn: \n\n "<< Inn << endl;
		debug << "\n A1inn: \n\n" << A1inn << endl;
		debug << "\n A1pinn: \n\n" << A1pinn << endl;	
		debug << "\n A1: \n\n" << A1 << endl;
		debug << "\n A1p: \n\n" << A1p << endl;
		debug << "\n A1top: \n\n" << A1top(1) << endl;
		debug << "\n A1ptop: \n\n" << A1ptop(1) << endl;
		debug << "\n tao: \n\n" << tao << endl;
		debug << "\n A:\n\n" << A << endl;		
	}			
}


//--------------------------------------------------------------------------------------------
 
void SAVE(void)								// Salva dati su file 

{
	if (imc <= 3)
	{
		ofstream inv_output(nomefile,ios::app);
		inv_output.setf(ios::fixed);
		inv_output.precision(4);
		inv_output.setf(ios::right);
		inv_output <<"\n";
		inv_output.width(6);
		inv_output << t;
		inv_output.width(30);
		inv_output << GDP1;
		inv_output.width(30);
		inv_output << GDP(1);
		inv_output.width(30);
		inv_output << EItot;
		inv_output.width(30);
		inv_output << SItot;
		//inv_output.width(20);
		//inv_output << Ntot;
		inv_output.width(25);
		inv_output << dNtot;
		inv_output.width(30);
		inv_output << Qtot1;
		inv_output.width(30);
		inv_output << Qtot2;
		inv_output.width(30);
		inv_output << Creal;
		inv_output.width(30);
		inv_output << Ir;
		inv_output.width(35);
		inv_output << Debtot1;
		inv_output.width(35);
		inv_output << Debtot2;
		inv_output.width(35);
		inv_output << Pitot1;
		inv_output.width(35);
		inv_output << Pitot2;
		inv_output.width(35);
		inv_output << Wtot1;
		inv_output.width(35);
		inv_output << Wtot2;
		inv_output.width(20);
		inv_output << LD;
		inv_output.width(15);
		inv_output << U(1);
		inv_output.width(45);
		inv_output << w(2);
		inv_output.width(15);
		inv_output << diff_w(2);
		//inv_output.width(15);
		//inv_output << ppi(1);
		inv_output.width(45);
		inv_output << cpi(1);
		inv_output.width(15);
		inv_output << diff_cpi(1);
		inv_output.width(20);
		inv_output << Am(1);
		inv_output.width(20);
		inv_output << A_sd;
		//inv_output.width(20);
		//inv_output << Am1;
		//inv_output.width(20);
		//inv_output << Am2;
		inv_output.width(15);
		inv_output << next1;
		inv_output.width(15);
		inv_output << next2;
		inv_output.width(25);
		inv_output << rw;
		inv_output.width(15);
		inv_output << Mutot;
		inv_output.width(50);
		inv_output << G(1);
		inv_output.width(50);
		inv_output << Tax;
		inv_output.width(50);
		inv_output << Def;
		inv_output.width(50);
		inv_output << Deb;
		inv_output.width(35);
		inv_output << DefonGDP;
		inv_output.width(35);
		inv_output << DebonGDP;
		inv_output.width(15);
		inv_output << H1;
		inv_output.width(15);
		inv_output << H2;
		inv_output.width(15);
		inv_output << N2lI;
		inv_output.width(15);
		inv_output << N2hI;
		inv_output.width(15);
		inv_output << N2lEI;
		inv_output.width(15);
		inv_output << N2hEI;
		inv_output.width(15);
		inv_output << DF1;
		inv_output.width(15);
		inv_output << DF2 << endl;
		inv_output.close();

 		ofstream inv_output1(nomefile1,ios::app);
 		inv_output1.setf(ios::fixed);
 		inv_output1.precision(4);
 		inv_output1.setf(ios::right);
 		inv_output1.width(6);
 		inv_output1 << t;
 		inv_output1.width(30);
 		inv_output1 << GDP1;
 		inv_output1.width(30);
 		inv_output1 << GDP(1);
 		inv_output1.width(30);
 		inv_output1 << EItot;
 		inv_output1.width(30);
 		inv_output1 << SItot;
 		inv_output1.width(25);
 		//inv_output1.width(20);
 		//inv_output1 << Ntot;
 		inv_output1 << dNtot;
 		inv_output1.width(30);
		inv_output1 << Qtot1;
		inv_output1.width(30);
		inv_output1 << Qtot2;
		inv_output1.width(30);
		inv_output1 << Creal;
		inv_output1.width(30);
		inv_output1 << Ir;
		inv_output1.width(35);
		inv_output1 << Debtot1;
		inv_output1.width(35);
		inv_output1 << Debtot2;
		inv_output1.width(35);
		inv_output1 << Pitot1;
		inv_output1.width(35);
		inv_output1 << Pitot2;
		inv_output1.width(35);
		inv_output1 << Wtot1;
		inv_output1.width(35);
		inv_output1 << Wtot2;
		inv_output1.width(20);
		inv_output1 << LD;
		inv_output1.width(15);
		inv_output1 << U(1);
		inv_output1.width(45);
		inv_output1 << w(2);
		inv_output1.width(15);
		inv_output1 << diff_w(2);
		//inv_output1.width(15);
		//inv_output1 << ppi(1);
		inv_output1.width(45);
		inv_output1 << cpi(1);
		inv_output1.width(15);
		inv_output1 << diff_cpi(1);
		inv_output1.width(20);
		inv_output1 << Am(1);
		inv_output1.width(20);
		inv_output1 << A_sd;
		//inv_output1.width(20);
		//inv_output1 << Am1;
		//inv_output1.width(20);
		//inv_output1 << Am2;
		inv_output1.width(15);
		inv_output1 << next1;
		inv_output1.width(15);
		inv_output1 << next2;
		inv_output1.width(25);
		inv_output1 << rw;
		inv_output1.width(15);
		inv_output1 << Mutot;
		inv_output1.width(50);
		inv_output1 << G(1);
		inv_output1.width(50);
		inv_output1 << Tax;
		inv_output1.width(50);
		inv_output1 << Def;
		inv_output1.width(50);
		inv_output1 << Deb;
		inv_output1.width(35);
		inv_output1 << DefonGDP;
		inv_output1.width(35);
		inv_output1 << DebonGDP;
		inv_output1.width(15);
		inv_output1 << H1;
		inv_output1.width(15);
		inv_output1 << H2;
		inv_output1.width(15);
		inv_output1 << N2lI;
		inv_output1.width(15);
		inv_output1 << N2hI;
		inv_output1.width(15);
		inv_output1 << N2lEI;
		inv_output1.width(15);
		inv_output1 << N2hEI;
		inv_output1.width(15);
		inv_output1 << DF1;
		inv_output1.width(15);
		inv_output1 << DF2 << endl;
 		inv_output1.close();
		}
    
	if (GDP(1) != Creal+Ir+dNtot && flagEXTRA!=1)
		cerr << "\n\n ERRORE contabilita nazionale!" << endl;
	
	ofstream inv_ymc(nomefile2,ios::app);
	inv_ymc.setf(ios::fixed);
	inv_ymc.precision(4);
	inv_ymc.setf(ios::right);
	inv_ymc.width(6);
	inv_ymc << t;
	inv_ymc.width(30);
	inv_ymc << GDP(1);
	inv_ymc.width(30);
	inv_ymc << Creal;
	inv_ymc.width(30);
	inv_ymc << Ir;
	inv_ymc.width(25);
	inv_ymc << dNtot;
	inv_ymc.width(30);
	inv_ymc << EItot;
	inv_ymc.width(30);
	inv_ymc << SItot;
	inv_ymc.width(25);
	inv_ymc << LD;
	inv_ymc.width(15);
	inv_ymc << U(1);
	inv_ymc.width(45);
	inv_ymc << w(2);
	inv_ymc.width(15);
	inv_ymc << diff_w(2);
	inv_ymc.width(45);
	inv_ymc << cpi(1);
	inv_ymc.width(15);
	inv_ymc << diff_cpi(1);
	inv_ymc.width(25);
	inv_ymc << Am(1);
	inv_ymc.width(25);
	inv_ymc << rw;
	inv_ymc.width(15);
	inv_ymc << Mutot;
	inv_ymc.width(50);
	inv_ymc << G(1);
	inv_ymc.width(50);
	inv_ymc << Tax;
	inv_ymc.width(50);
	inv_ymc << Def;
	inv_ymc.width(50);
	inv_ymc << Deb;
	inv_ymc.width(35);
	inv_ymc << DefonGDP;
	inv_ymc.width(35);
	inv_ymc << DebonGDP;
	inv_ymc.width(15);
	inv_ymc << H1;
	inv_ymc.width(15);
	inv_ymc << H2;
	inv_ymc.width(15);
	inv_ymc << DF1;
	inv_ymc.width(15);
	inv_ymc << DF2 << endl;
	inv_ymc.close();
}

//--------------------------------------------------------------------------------------------

void UPDATE(void)							// Aggiorno variabili per prossimo turno
 {
	 //cout << "\nUpdate..." <<endl;	 
	 Em2(2)=Em2(1);
	 Am(2)=Am(1);
	 cpi(2)=cpi(1);
	 ppi(2)=ppi(1);
	 diff_cpi(2)=diff_cpi(1);
	 U(2)=U(1);
	 w(3)=w(2);
	 w(2)=w(1);
	 wstar(2)=wstar(1);
	 diff_w(2)=diff_w(1);
	 GDP(3)=GDP(2);
	 GDP(2)=GDP(1);
	 G(2)=G(1);
	 Divtot(2)=Divtot(1);
	 Ir2=Ir;
	 dNtot2=dNtot;
	 A1front(2)=A1front(1);
	 A1pfront(2)=A1pfront(1);
	 A1top(2)=A1top(1);
	 A1ptop(2)=A1ptop(1);
	 
	 for (i=1; i<=N1; i++)
	 {
		 W1(2,i)=W1(1,i);
		 S1(2,i)=S1(1,i);
		 Deb1(2,i)=Deb1(1,i);
		 p1(2,i)=p1(1,i);
		 RD(2,i)=RD(1,i);
		 size1(2,i)=size1(1,i);
		 f1(2,i)=f1(1,i);
	 }

	 for (j=1; j<=N2; j++)
	 {
		 E2(2,j)=E2(1,j);
		 f2(3,j)=f2(2,j);
		 f2(2,j)=f2(1,j);
		 W2(2,j)=W2(1,j);
		 D2(5,j)=D2(4,j);
		 D2(4,j)=D2(3,j);
		 D2(3,j)=D2(2,j);
		 D2(2,j)=D2(1,j);
		 De(3,j)=De(2,j);
		 De(2,j)=De(1,j);
		 Q2(2,j)=Q2(1,j);
		 c2(2,j)=c2(1,j);
		 N(2,j)=N(1,j);
		 Nm(2,j)=Nm(1,j);
		 I(2,j)=I(1,j);
		 EI(2,j)=EI(1,j);
		 SI(2,j)=SI(1,j);
		 S2(2,j)=S2(1,j);
		 Deb2(2,j)=Deb2(1,j);
		 p2(2,j)=p2(1,j);
		 mu2(2,j)=mu2(1,j);
		 size2(2,j)=size2(1,j);

		 for (i=1; i<=N1; i++)						// Aumento et‚àö‚Ä† macchine
		 {												
		 	for (tt=t0; tt<=t; tt++)
			{
				if (gtemp[tt-1][i-1][j-1] > 0)
					age[tt-1][i-1][j-1]++;
			}
		 }
	 }
 }
//--------------------------------------------------------------------------------------------
 
void OVERBOOST(void)

 {
	 t00=t0;
	 flag=0;
	 mintao=tao.Minimum();
	 for (tt=t00; tt<=t && flag == 0; tt++)
	 {	 
		 for (i=1; i<=N1 && flag == 0; i++)
		 { 
			 for (j=1; j<=N2 && flag == 0; j++)
			 {
				 if (g[tt-1][i-1][j-1] > 0 || gtemp[tt-1][i-1][j-1] > 0)
					 flag=1;
			 }
		 }
		 if (flag == 0)
			 t0++;
	 }
 }

//--------------------------------------------------------------------------------------------

 double ROUND(double x)						// Arrotondamento
 {
	 double x_floor=floor(x);
	 double resto=x-x_floor;
	 if (resto > 0.5) x=x_floor+1;
	 else x=x_floor;
	 return x;
 }

//--------------------------------------------------------------------------------------------

string LONGSTR(long n)						// Convert int to string
{
	stringstream s;
	s << n;									//scrive su string stream
	return s.str();
}

//--------------------------------------------------------------------------------------------

void GENFILEOUTPUT(char s[],long number)				    // Writes on FileName[] the string: 
													// s + number + ".txt"
{
	strcpy(nomefile,s);								// Copia s in nomefile
	string ss1 = LONGSTR(number);					// Converte numero in stringa 	
	const char* p1 = ss1.c_str();					
	char* name=strcat(nomefile,p1);					// Concatena 2 stringhe
	strcat(nomefile,".txt");
}

//-------------------------------------------------------------------------------------------

void GENFILEOUTPUT1(char s1[],long number1)				    // Writes on FileName[] the string: 
													// s + number + ".txt"
{
	strcpy(nomefile1,s1);								// Copia s in nomefile
	string ss2 = LONGSTR(number1);					// Converte numero in stringa 	
	const char* p2 = ss2.c_str();					
	char* name1=strcat(nomefile1,p2);					// Concatena 2 stringhe
	strcat(nomefile1,".txt");
}

//-------------------------------------------------------------------------------------------

void GENFILESIZE1(char s0[],long number0)				    // Writes on FileName[] the string: 
															// s + number + ".txt"
{
	strcpy(nomefile0,s0);								// Copia s in nomefile
	string ss0 = LONGSTR(number0);					// Converte numero in stringa 	
	const char* p0 = ss0.c_str();					
	char* name0=strcat(nomefile0,p0);					// Concatena 2 stringhe
	strcat(nomefile0,".txt");
}
//-------------------------------------------------------------------------------------------

void GENFILESIZE2(char s3[],long number3)				    // Writes on FileName[] the string: 
													// s + number + ".txt"
{
	strcpy(nomefile3,s3);								// Copia s in nomefile
	string ss4 = LONGSTR(number3);					// Converte numero in stringa 	
	const char* p4 = ss4.c_str();					
	char* name3=strcat(nomefile3,p4);					// Concatena 2 stringhe
	strcat(nomefile3,".txt");
}

//--------------------------------------------------------------------------------------------

void GENFILEYMC(char s2[],long number2)				    // Writes on FileName[] the string: 
													// s + number + ".txt"
{
	strcpy(nomefile2,s2);								// Copia s in nomefile
	string ss2 = LONGSTR(number2);					// Converte numero in stringa 	
	const char* p3 = ss2.c_str();					
	char* name2=strcat(nomefile2,p3);					// Concatena 2 stringhe
	strcat(nomefile2,".txt");
}

//-----------------------------------------------------------------------------------------

void GENFILEPROD2(char s4[],long number4)				    // Writes on FileName[] the string: 
													// s + number + ".txt"
{
	strcpy(nomefile4,s4);								// Copia s in nomefile
	string ss5 = LONGSTR(number4);					// Converte numero in stringa 	
	const char* p5 = ss5.c_str();					
	char* name4=strcat(nomefile4,p5);					// Concatena 2 stringhe
	strcat(nomefile4,".txt");
}

//-----------------------------------------------------------------------------------------

void GENFILEPROD1(char s5[],long number5)				    // Writes on FileName[] the string: 
															// s + number + ".txt"
{
	strcpy(nomefile5,s5);								// Copia s in nomefile
	string ss6 = LONGSTR(number5);					// Converte numero in stringa 	
	const char* p6 = ss6.c_str();					
	char* name5=strcat(nomefile5,p6);					// Concatena 2 stringhe
	strcat(nomefile5,".txt");
}

//-----------------------------------------------------------------------------------------

void GENFILEPRODALL1(char s6[],long number6)				    // Writes on FileName[] the string: 
															// s + number + ".txt"
{
	strcpy(nomefile6,s6);								// Copia s in nomefile
	string ss7 = LONGSTR(number6);					// Converte numero in stringa 	
	const char* p7 = ss7.c_str();					
	char* name6=strcat(nomefile6,p7);					// Concatena 2 stringhe
	strcat(nomefile6,".txt");
}

//-----------------------------------------------------------------------------------------

void GENFILEPRODALL2(char s7[],long number7)				    // Writes on FileName[] the string: 
															// s + number + ".txt"
{
	strcpy(nomefile7,s7);								// Copia s in nomefile
	string ss8 = LONGSTR(number7);					// Converte numero in stringa 	
	const char* p8 = ss8.c_str();					
	char* name7=strcat(nomefile7,p8);					// Concatena 2 stringhe
	strcat(nomefile7,".txt");
}
